Fractional Boundary Element Analysis of Nonlinear Three-Temperature Thermoelectric Problems

1. Introduction

Utilizing initial conditions and boundary conditions at m and Max Planck proposed the quantum theory of electromagnetic radiation in 1900, and Albert Einstein proposed the concept of photons and phonons in 1905 and 1907 to explain why temperature varies with the specific heat of solid crystals. A phonon is a quantum of vibrational mechanical energy produced by a lattice of oscillating atoms, where the thermal energy of the atoms causes the lattice to vibrate. This generates compression mechanical waves, which carry heat and sound through the anisotropic material. The careful study of phonons is an important part of solid-state physics because it is important for many thermal and acoustic properties of solids, as well as a variety of low- and high-temperature superconductive technologies.

Many bodies undergo significant changes in their properties as a result of the application of an electric or magnetic field, allowing for the profitable application of this fact in technological applications. The development of electrorheological and magnetic fluids has piqued the interest of many people due to their potential applications in clutches, actuators, shock absorbers, valves, and exercise equipment, to name a few. Field dependent solids, also known as electro-active elastomers, were recently developed by infusing elastomers with electrorheological fluids or embedding them with electrically conducting particles. Magneto-active elastomers, on the other hand, were created by encapsulating elastomers with magnetically responsive particles. Such field dependent solids have potential applications due to the change in structure and the resulting effect on the elasticity and compliance of the material, such field dependent solids have potential applications in a wide range of applications.

Adaptive algorithms based on fractional calculus (FC) have been proposed for parameter estimation of various problems [1,2], and they outperform standard adaptive methodologies in terms of convergence speed and estimation accuracy [3]. FC is concerned with the calculation of real-order derivatives and integrals [4]. The FC has been successfully used to solve problems in circuit design [5], artistic paintings [6], vibration analysis [7], hydro turbine systems [8], control engineering [9,10], nanotechnology [11], and biological processes [12,13]. Fractional adaptive algorithms are used to estimate the parameters of power signals. The fractional order least mean square (FOLMS) technique is used to estimate the amplitude and phase of power signals [14]. The momentum FOLMS (mFOLMS) was created to speed up the convergence speed of conventional FOLMS [15]. Experiments have shown that the mFOLMS algorithm outperforms the traditional LMS and standard FOLMS algorithms in terms of convergence speed. The FOLMS and mFOLMS performance are unaffected by the fractional order. The so-called innovative FOLMS (I-FOLMS) adaptive algorithm [16], on the other hand, was recently developed, indicating a strong reliance on fractional order. The I-FOLMS solves the opposing requirements of fast convergence speed and small steady-state error by using an appropriate fractional order selection. The fractional order parameter was discovered to be capable of influencing the initial convergence speed and estimation accuracy of I-FOLMS. These intriguing properties of the I-FOLMS necessitate further research into the I-FOLMS for parameter estimation of power signals, as well as comparisons with the regular FOLMS.

Knopoff [17] and Chadwick [18], followed by Kaliski and Petykiewicz [19], pioneered the magnetoelasticity foundations. Magnetic theory development and application were previously solely based on magnetic experiments. Because of the rapid development of high-performance computational methods and computer hardware, efficient and accurate computational methods for modeling and simulation of real magnetic experiments have been used, particularly when the magnetic experiment is difficult, dangerous, or expensive.

The interaction between magnetic and strain fields in a thermoelastic solid is receiving increased attention due to its numerous applications in geophysics, plasma physics, and other fields. All of the articles cited above assumed that interactions between the two fields occur via Lorentz forces appearing in equations of motion and a term entering Ohm’s law and representing the electric field created by the velocity of a material particle traveling in a magnetic field. In these investigations, the heat equation under consideration is typically the uncoupled or coupled theory, rather than the generalised one. Ezzat and Awad developed a model of micropolar generalised magneto-thermoelasticity based on modified Ohm’s and Fourier’s laws [20]. In the literature, several fractional-order models have been investigated for various applications. There is no general analytical solution due to the computational difficulty of solving complex fractional nonlinear thermoelasticity problems. Numerical methods, including the BEM, should be used to solve such problems. Different three-temperature theories have been investigated in the context of micropolar-thermoelasticity [21], carbon nanotube fiber reinforced composites [22], micropolar piezothermoelasticity [23], Micropolar Magneto-thermoviscoelasticity [24] and Magneto-thermooviscoelasticity [25].

In this paper, we presented a new fractional order theory of functionally graded magnetic thermoelectric materials. This theory’s application to three-temperature nonlinear generalized thermoelasticity is solved using boundary element analysis. Because of its advantages, such as dealing with more complex shapes of functionally graded magnetic thermoelectric materials and not requiring the discretization of the internal domain, BEM has low RAM and CPU usage. As a result, it is adaptable and effective for dealing with complex nonlinear 3T thermoelectric FGM problems. The numerical results demonstrate the effects of a magnetic field and a graded parameter on thermal stresses in FGMTM. The numerical results also confirm the validity and accuracy of the proposed technique.

2. Formulation of the Problem

A Cartesian coordinate system for a 2D functionally graded magnetic thermoelectric structure with thickness of h and length $L$ as shown in Figure 1. It is assumed to be subjected to an electric potential Φ(x, z, t) along the $0\mathrm{x}$ direction and placed in an external constant magnetic field ${\mathrm{H}}_0$ within the region $\mathrm{R}=\left\{0<\mathrm{x}<L,0<\mathrm{y}<h\right\}$ which bounded by boundary $\mathrm{S}$, where ${\mathrm{S}}_{\mathrm{i}}$ $\left(\mathrm{i}=1,2,3,4\right)$ are subsets of $\mathrm{S}$ such that ${\mathrm{S}}_1+{\mathrm{S}}_2={\mathrm{S}}_3+{\mathrm{S}}_4=\mathrm{S}$.

The governing equations for fractional order three-temperature nonlinear generalized thermoelastic problems of functionally graded magnetic thermoelectric materials can be written as [26]

$${\sigma }_{ij,j}+{\mu }_0{\left(x+1\right)}^m{\epsilon }_{ijk} J_k{H}_j=\rho {\left(x+1\right)}^m{\ddot{u}}_i$$

(1)

where

$${\sigma }_{ij}={\left(x+1\right)}^m\left[C_{ijkl} e{\delta }_{ij}+\check{\alpha } \left(u_{j,i}-{\epsilon }_{ijk}{\omega }_k\right)-{\beta }_{ij} \left(\theta +{\tau }_1\dot{\theta }\right)\right]$$

(2)

$$J_i={\sigma }_0\left(E_i+{\epsilon }_{ijk}{\dot{u}}_k{B}_j-k_0{T}_{,i}\right)$$

(3)

$$B_i={\mu }_0{H}_i$$

(4)

The fractional order three-temperature radiative heat conduction equations coupled with electron, ion and phonon temperatures can be expressed as

$$D_{\tau }^a{T}_{\alpha }\left(r, \tau \right)=\xi \nabla \left[{\mathbb{K}}_{\alpha } \nabla T_{\alpha }(r, \tau )\right]+\xi \overline{\mathbb{W}}\left(r, \tau \right), \xi ={\displaystyle \frac{1}{c_{\alpha }\rho {\delta }_1}}$$

(5)

where

$$\overline{\mathbb{W}}\left(r, \tau \right)=\{\begin{array}{c}-\rho {\mathbb{W}}_{ei} \left(T_e-T_i\right)-{\rho \mathbb{W}}_{ep} \left(T_e-T_p\right)+\stackrel{̿}{\mathbb{W}}\left(r, \tau \right), \alpha =e, {\delta }_1=1 \\ \rho {\mathbb{W}}_{ei} \left(T_e-T_i\right)+\stackrel{̿}{\mathbb{W}}\left(r, \tau \right), \alpha =i,{ \delta }_1=1 \\ \rho {\mathbb{W}}_{ep} \left(T_e-T_p\right)+\stackrel{̿}{\mathbb{W}}\left(r, \tau \right), \alpha =p, {\delta }_1={\displaystyle \frac{4}{\rho }}{T}_p^3 \end{array}$$

(6)

in which

$$\stackrel{̿}{\mathbb{W}}\left(r, \tau \right)=-{\delta }_{2j}{\mathbb{K}}_{\alpha } {\dot{T}}_{\alpha ,ab}+{\beta }_{ab}{T}_{\alpha 0}\left[Å{\delta }_{1j}{\dot{u}}_{a,b}+\left({\tau }_0+{\delta }_{2j}\right){\ddot{u}}_{a,b}\right]$$

$$+\rho c_{\alpha }\left[\left({\tau }_0+{\delta }_{1j}{\tau }_2+{\delta }_{2j}\right){\ddot{T}}_{\alpha }\right]-\rho {\pi }_0{ J}_{i,j}$$

And

$${\mathbb{W}}_{ei}=\rho {\mathbb{A}}_{ei}{T}_e^{-2/3}, {\mathbb{W}}_{ep}=\rho {\mathbb{A}}_{ep}{T}_e^{-1/2}, {\mathbb{K}}_{\alpha }={\mathbb{A}}_{\alpha }{T}_{\alpha }^{5/2}, \alpha =e, i, {\mathbb{K}}_p={\mathbb{A}}_p{T}_p^{3+\mathbb{B}}$$

(7)

The total energy of a unit mass can be described as follows:

$$P=P_e+P_i+P_p, P_e=c_e{T}_e, P_i=c_i{T}_i,P_p={\displaystyle \frac{1}{\rho }}{c}_p{T}_p^4$$

(8)

Initial and boundary conditions can be written as

$$T_{\alpha }\left(x, y, 0\right)=T_{\alpha }^0\left(x, y\right)=g_1\left(x, \tau \right)$$

(9a)

$${{\mathbb{K}}_{\alpha }{\displaystyle \frac{\partial T_{\alpha }}{\partial n}}|}_{{\mathsf{\Gamma }}_1}=0, \alpha =e, i, {T_p|}_{{\mathsf{\Gamma }}_1}=g_2\left(x, \tau \right)$$

(9b)

$${{\mathbb{K}}_{\alpha }{\displaystyle \frac{\partial T_{\alpha }}{\partial n}}|}_{{\mathsf{\Gamma }}_2}=0, \alpha =e, i, p$$

(9c)

3. BEM Simulation for Temperature Field

The boundary element method is used in this section to solve the nonlinear time-dependent two dimensions three temperature (2D-3T) radiation diffusion equations that are coupled by electron, ion, and photon temperatures.

Caputo’s integral definition and Grunwald-Letnikov integral definition are consistent with the Riemann-Liouville integral definition. Also, Caputo’s derivative definition and Grunwald-Letnikov derivative definition are consistent with the Riemann-Liouville derivative definition. When the Riemann-Liouville or Grunwald-Letnikov definitions are compared to the Caputo definition, the functions that are derivable in the Caputo sense are significantly fewer. According to finite difference scheme of Caputo at times $\left(f+1\right)∆\tau$ and $f∆\tau$, we obtain [27]

$$D_{\tau }^a{T}_{\alpha }^{f+1}+D_{\tau }^a{T}_{\alpha }^f\approx {\displaystyle \sum _{j=0}^k}{W}_{a,j}\left(T_{\alpha }^{f+1-j}\left(r\right)-T_{\alpha }^{f-j}\left(r\right)\right),\left(f=1, 2, \dots ., F\right)$$

(10)

where

$$W_{a,0}={\displaystyle \frac{{\left(∆\tau \right)}^{-a}}{\Gamma \left(2-a\right)}} , W_{a,j}=W_{a,0}\left({\left(j+1\right)}^{1-a}-{\left(j-1\right)}^{1-a}\right), j=1, 2, \dots ., F$$

(11)

Based on equation (10), the fractional order heat equations (5) can be replaced by the following system

$$W_{a,0}{T}_{\alpha }^{f+1}\left(r\right)-{\mathbb{K}}_{\alpha }\left(x\right)T_{\alpha ,II}^{f+1}\left(r\right)-{{\mathbb{K}}_{\alpha }}_{,I}\left(x\right)T_{\alpha ,I}^{f+1}\left(r\right)$$

$$=W_{a,0}{T}_{\alpha }^f\left(r\right)-{\mathbb{K}}_{\alpha }\left(x\right)T_{\alpha ,II}^f \left(r\right)-{{\mathbb{K}}_{\alpha }}_{,I}\left(x\right)T_{\alpha ,J}^f \left(r\right)$$

$$-{\displaystyle \sum _{j=1}^f}{W}_{a,j}\left(T_{\alpha }^{f+1-j}\left(r\right)-T_{\alpha }^{f-j}\left(r\right)\right)+{\overline{\mathbb{W}}}_m^{f+1}\left(x, \tau \right)$$

$$+{\overline{\mathbb{W}}}_m^f\left(x, \tau \right), f=0, 1, 2, \dots , F$$

(12)

Based on the fundamental solution of (12), the direct formulation of boundary integral equation corresponding to (5) can be expressed as

$$C{T}_{\alpha }={\displaystyle \frac{D}{{\mathbb{K}}_{\alpha }}}{\int }_O^{\tau }{\int }_S\left[T_{\alpha }{q}^{* }-T_{\alpha }^{*}q\right]dS d\tau +{\displaystyle \frac{D}{{\mathbb{K}}_{\alpha }}}{\int }_O^{\tau }{\int }_{R}b{T}_{\alpha }^{*}dR d\tau +{\int }_R{T_{\alpha }^i{T}_{\alpha }^{*}|}_{\tau =0}dR$$

(13)

which can be written in the absence of internal heat sources as follows

$$C{T}_{\alpha }={\int }_S\left[T_{\alpha }{q}^{*}-T_{\alpha }^{*}q\right] dS-{\int }_R{\displaystyle \frac{{\mathbb{K}}_{\alpha }}{D}} {\displaystyle \frac{\partial T_{\alpha }^{*}}{\partial \tau }}{T}_{\alpha } dR$$

(14)

We assume that the time derivative of temperature can be approximated by a series of known functions in order to transform the domain integral in (14) to the boundary. $f^j\left(r\right)$ and unknown coefficients $a^j\left(\tau \right)$ as

$${\displaystyle \frac{\partial T_{\alpha }}{\partial \tau }}\cong {\displaystyle \sum _{j=1}^N}{f}^j{\left(r\right)}^j{a}^j\left(\tau \right)$$

(15)

Also, we assume that ${\widehat{T}}_{\alpha }^j$ is a solution of

$${\nabla }^2{\widehat{T}}_{\alpha }^j=f^j$$

(16)

Thus, equation (14) results in the following boundary integral equation

$$C T={\int }_S\left[T_{\alpha }{q}^{*}-T_{\alpha }^{*}q\right] dS+{\displaystyle \sum _{j=1}^N}{a}^j\left(\tau \right)D^{-1}\left(C{\widehat{T}}_{\alpha }^j-{\int }_S\left[T_{\alpha }^j{q}^{*}-{\widehat{q}}^j{T}_{\alpha }^{*}\right] dS\right)$$

(17)

where

$${\widehat{q}}^j=-{\mathbb{K}}_{\alpha }{\displaystyle \frac{\partial {\widehat{T}}_{\alpha }^j}{\partial n}}$$

(18)

and

$$a^j\left(\tau \right)={\displaystyle \sum _{i=1}^N}{f}_{ji}^{-1}{\displaystyle \frac{\partial T\left(r_i, \tau \right)}{\partial \tau }}$$

(19)

In which the entries of $f_{ji}^{-1}$ are the coefficients of $F^{-1}$ with matrix $F$ defined as [24]

$${\left\{F\right\}}_{ji}=f^j\left(r_i\right)$$

(20)

Using the standard boundary element discretization scheme for equation (17) and using equation (19), we obtain the following set of ordinary differential equations

$$C {\dot{T}}_{\alpha }+H T_{\alpha }=G Q$$

(21)

where matrices $H$ and $G$ are depending on the current time step, boundary geometry and material properties, $T_{\alpha }$ and $Q$ are, respectively, temperature and heat flux vectors at boundary nodes, and b is the internal heat generation vector

The diffusion matrix can be defined as

$$C=-\left[H {\widehat{T}}_{\alpha }-G \widehat{Q}\right]F^{-1}{D}^{-1}$$

(22)

with

$${\left\{\widehat{T}\right\}}_{ij}={\widehat{T}}^j\left(x_i\right)$$

(23)

$${\left\{\widehat{Q}\right\}}_{ij}={\widehat{q}}^j\left(x_i\right)$$

(24)

In order to solve equation (21) numerically the functions $T_{\alpha }$ and $q$ are interpolated as

$$T_{\alpha }=\left(1-\theta \right) T_{\alpha }^m+\theta T_{\alpha }^{m+1}$$

(25)

$$q=\left(1-\theta \right)q^m+\theta q^{m+1}$$

(26)

where the parameter $\theta ={\displaystyle \frac{\tau -{\tau }^m}{{\tau }^{m+1}-{\tau }^m}} , 0\le \theta \le 1$ determines the practical time $\tau$ in the present time step.

By differentiating equation (25) with respect to time we get

$${\dot{T}}_{\alpha }={\displaystyle \frac{d{T}_{\alpha }}{d\theta }}{\displaystyle \frac{d\theta }{d\tau }}={\displaystyle \frac{T_{\alpha }^{m+1}-T_{\alpha }^m}{{\tau }^{m+1}-{\tau }^m}}={\displaystyle \frac{T_{\alpha }^{m+1}-T_{\alpha }^m}{∆{\tau }^m}}$$

(27)

By substituting from Equations (25) - (27) into Equation (21), we obtain

$$\left({\displaystyle \frac{C}{∆{\tau }^m}}+\theta H\right)T_{\alpha }^{m+1}-\theta G{Q}^{m+1}=\left({\displaystyle \frac{C}{∆{\tau }^m}}-\left(1-\theta \right)H\right)T_{\alpha }^m+\left(1-\theta \right)G{Q}^m$$

(28)

By using initial and boundary conditions at $∆{\tau }^m$ and considering the previous time step solution as initial values for the next step, we obtain the following linear algebraic system

$$\mathbb{a}\mathsf{{\rm X}}=\mathbb{b}$$

(29)

where $\mathbb{a}$ is unknown matrix, and $\mathsf{{\rm X}}$ and $\mathbb{b}$ are known matrices

4. BEM Simulation for Displacement and Microrotation Fields

Using the weighted residual method, the governing equations (1) and (2) can be transformed into the following integral equations

$${\int }_R\left({\sigma }_{ij,j}+U_i\right) u_i^{*} dR=0$$

(30)

in which

$$U_i={\mu }_0{\left(x+1\right)}^m{\epsilon }_{ijk}{J}_k{H}_j-\rho {\ddot{u}}_i$$

(31)

The boundary conditions are

$$u_i={\overline{u}}_i \mathrm{o}\mathrm{n} S_1$$

(32)

$${\lambda }_i={\sigma }_{ij}{n}_j={\overline{\lambda }}_i \mathrm{o}\mathrm{n} S_2$$

(33)

The integration of the first term of equations (33) and (34) leads to

$$-{\int }_R{\sigma }_{ij} u_{i,j}^{*} dR+{\int }_R{U}_i u_i^{*} dR=-{\int }_{S_2}{\lambda }_i u_i^{*} dS$$

(34)

According to Huang and Liang [28], the boundary integral equation can be written as

$$-{\int }_R{\sigma }_{ij,j} u_i^{*} dR+{\int }_R{U}_i u_i^{*} dR={\int }_{S_2} \left({\lambda }_i-{\overline{\lambda }}_i\right)u_i^{*} dS+{\int }_{S_1}\left({\overline{u}}_i-u_i\right) {\lambda }_i^{*} dS$$

(35)

The integration of (35)’s left-hand side by parts results in

$$-{\int }_R{\sigma }_{ij} {\epsilon }_{ij}^{*} dR+{\int }_R{U}_i u_i^{*} dR=-{\int }_{S_2} {\overline{\lambda }}_i{ u}_i^{*} dS-{\int }_{S_1} {\lambda }_i u_i^{*} dS+{\int }_{S_1}\left({\overline{u}}_i-u_i\right) {\lambda }_i^{*} dS$$

(36)

According to Eringen [29], the elastic stress can be written as

$${\sigma }_{ij}={\mathbb{A}}_{ijkl} {\epsilon }_{kl}, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} {\mathbb{A}}_{ijkl}={\mathbb{A}}_{klij}$$

(37)

Hence, equation (36) may be expressed as

$$-{\int }_R{\sigma }_{ij}^{*} {\epsilon }_{ij} dR+{\int }_R{U}_i u_i^{*} dR=-{\int }_{S_2} {\overline{\lambda }}_i{ u}_i^{*} dS-{\int }_{S_1} {\lambda }_i u_i^{*} dS+{\int }_{S_1}\left({\overline{u}}_i-u_i\right) {\lambda }_i^{*} dS$$

(38)

By applying the integration by parts again to the left-hand side of (38), we obtain

$${\int }_R {\sigma }_{ij,j}^{*} u_i dR=-{\int }_S{u}_i^{*} {\lambda }_i dS+{\int }_S{\lambda }_i^{*} u_i dS$$

(39)

where

$${\sigma }_{lj,j}^{*}+{∆}^n{e}_l=0$$

(40)

According to Dragos [30], the fundamental solution may be written as

$$u_i^{*}=u_{li}^{*} e_l,{\lambda }_i^{*}={\lambda }_{li}^{*} e_l,$$

(41)

The weighting functions for $U_i=0$ and $V_i={∆}^n$ may be expressed as follows:

$${\sigma }_{ij,j}^{**}=0$$

(42)

On the basis of Dragos [30], the fundamental solution may be expressed as

$$u_i^{*}=u_{li}^{**} e_l,{\lambda }_i^{*}={\lambda }_{li}^{**} e_l$$

(43)

Using the above two sets of weighting functions into (39) we have

$$C_{li}^n{u}_i^n=-{\int }_S{\lambda }_{li}^{*} u_{i}dS+{\int }_S{u}_{li}^{*} {\lambda }_i dS$$

(44)

$$C_{li}^n{\omega }_i^n=-{\int }_S{\lambda }_{li}^{**} u_{i}dS+{\int }_S{u}_{li}^{**} {\lambda }_i dS$$

(45)

Thus, we obtain

$$C^n{\mathbb{q}}^n=-{\int }_S{\mathbb{p}}^{*}\mathbb{q}dS+{\int }_S{\mathbb{q}}^{*}\mathbb{p}dS$$

(46)

in which

$$C^n=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right], {\mathbb{q}}^{*}=\left[\begin{array}{cc}{u}_{11}^{*}& u_{12}^{*}\\ u_{21}^{*}& u_{22}^{*}\end{array}\right], {\mathbb{p}}^{*}=\left[\begin{array}{cc}{\lambda }_{11}^{*}& {\lambda }_{12}^{*}\\ {\lambda }_{21}^{*}& {\lambda }_{22}^{*}\end{array}\right],\mathbb{q}=\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right],\mathbb{p}=\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\end{array}\right]$$

(47)

Now, we introduce the following relations

$$\mathbb{q}=\psi {\mathbb{q}}^j,\mathbb{p}=\psi {\mathbb{p}}^j$$

(48)

By discretizing the boundary, we can write (46) as

$$C^n{\mathbb{q}}^n={\displaystyle \sum _{j=1}^{N_e}}\left[-{\int }_{{\mathsf{\Gamma }}_j}{\mathbb{p}}^{*}\psi d\mathsf{\Gamma }\right]{\mathbb{q}}^j+{\displaystyle \sum _{j=1}^{N_e}}\left[{\int }_{{\mathsf{\Gamma }}_j}{\mathbb{q}}^{*}\psi d\mathsf{\Gamma }\right]{\mathbb{p}}^j$$

(49)

Which can be expressed as

$$C^i{\mathbb{q}}^i=-{\displaystyle \sum _{j=1}^{N_e}}{\widehat{\mathbb{H}}}^{ij}{\mathbb{q}}^j+{\displaystyle \sum _{j=1}^{N_e}}{\widehat{\mathbb{G}}}^{ij}{\mathbb{p}}^j$$

(50)

By employing the following formula

$${\mathbb{H}}^{ij}=\{\begin{array}{c}{\widehat{\mathbb{H}}}^{ij} if i\ne j\\ {\widehat{\mathbb{H}}}^{ij}+C^i if i=j\end{array}$$

(51)

Hence, equation (50) may be expressed as

$${\displaystyle \sum _{j=1}^{N_e}}{\mathbb{H}}^{ij}{\mathbb{q}}^j={\displaystyle \sum _{j=1}^{N_e}}{\widehat{\mathbb{G}}}^{ij}{\mathbb{p}}^j$$

(52)

The global matrix system equation for all $i$ nodes can be written as follows

$$\mathbb{H}\mathbb{Q}=\mathbb{G}\mathbb{P}$$

(53)

where $\mathbb{Q}$ denotes the displacements and $\mathbb{P}$ denotes the tractions.

Now, we can write (53) into the following form

$$\mathbb{A}\mathbb{X}=\mathbb{B}$$

(54)

In Matlab (R2018a), an explicit staggered predictor-corrector scheme based on the communication-avoiding generalized minimal residual (CA-GMRES) method is efficiently implemented for solving the resulting simultaneous linear algebraic systems to obtain the temperature and displacement fields.

5. Numerical Results and Discussion

In the context of functionally graded magnetic thermoelectric materials, the proposed BEM technique can be applied to a wide range of fractional-order nonlinear generalised thermal stress problems. The BEM discretization was performed using 42 boundary elements and 68 internal points, as shown in Figure 2.

Figure 3. Shows the variation of the thermal stress ${\sigma }_{11}$ along $x$-axis for different values of fractional order parameter ($a=0.3, 0.6 and 0.9$). It is shown from this figure that the thermal stress ${\sigma }_{11}$ decreases with the increase of $x$ until $x=0.9$. Then it increases with the increase of $x$. It is also shown from this figure that the thermal stress ${\sigma }_{11}$ decreases with the increase of fractional order parameter.

Figure 4. Shows the variation of the thermal stress ${\sigma }_{12}$ along $x$-axis for different values of fractional order parameter ($a=0.3, 0.6 and 0.9$). It is shown from this figure that the thermal stress ${\sigma }_{12}$ decreases with the increase of $x$. It is also shown from this figure that the thermal stress ${\sigma }_{21}$ increases with the increase of fractional order parameter.

Figure 5 shows the variation of the thermal stress ${\sigma }_{22}$ along $x$-axis for different values of fractional order parameter ($a=0.3, 0.6 and 0.9$). It is shown from this figure that the thermal stress ${\sigma }_{22}$ increases with the increase of $x$. It is also shown from this figure that the thermal stress ${\sigma }_{22}$ decreases with the increase of fractional order parameter.

Figure 6. Shows the variation of the thermal stress ${\sigma }_{11}$ along $x$-axis for different values of functionally graded parameter ($m=0.3, 0.6 and 0.9$). It is shown from this figure that the thermal stress ${\sigma }_{11}$ increases with the increase of $x$ until $x=0.9$. Then it decreases with the increase of $x$. It is also shown from this figure that the thermal stress ${\sigma }_{11}$ increases with the increase of functionally graded parameter.

Figure 7. Shows the variation of the thermal stress ${\sigma }_{12}$ along $x$-axis for different values of functionally graded parameter ($m=0.3, 0.6 and 0.9$). It is shown from this figure that the thermal stress ${\sigma }_{12}$ increases with the increase of $x$. It is also shown from this figure that the thermal stress ${\sigma }_{21}$ decreases with the increase of functionally graded parameter.

Figure 8 Shows the variation of the thermal stress ${\sigma }_{22}$ along $x$-axis for different values of functionally graded parameter ($m=0.3, 0.6 and 0.9$). It is shown from this figure that the thermal stress ${\sigma }_{22}$ increases with the increase of $x$. It is also shown from this figure that the thermal stress ${\sigma }_{22}$ decreases with the increase of functionally graded parameter.

It is noted from figures 3-8 that the fractional order parameter and functionally graded parameter have a strong effect on the thermal stress ${\mathsf{\sigma }}_{11}$, ${\mathsf{\sigma }}_{12}$ and ${\mathsf{\sigma }}_{22}$ in the functionally graded magnetic thermoelectric materials.

Table 1 shows a comparison of required computer resources for the current BEM results, GFDM results of Gu et al. [31] and FEM–NMM results of An et al. [32] of modeling of fractional nonlinear three-temperature (3T) thermoelectric problems.

There were no published results to demonstrate the validity of the proposed technique’s results. Some literatures, on the other hand, can be regarded as special cases of the considered general study. Figure 9 shows the variation of the special case thermal stress ${\sigma }_{11}$ along $x$-axis for BEM, generalized finite difference method (GFDM) and combined finite element method/normal mode method (FEM–NMM) in the case of fractional order $\left(a=0.6\right)$ homogeneous $\left(m=0.0\right)$. Figure 10 shows the propagation of the thermal stress ${\sigma }_{11}$ along $x$-axis for BEM, GFDM and FEM–NMM in the case of zero fractional order, $\left(a=0.0\right)$ and non-homogeneous $\left(m=0.6\right)$. These findings for thermal stress ${\sigma }_{11}$ in functionally graded magnetic thermoelectric materials, show that the BEM findings are in a very good agreement with the GFDM findings of Gu et al. [31], and FEM–NMM findings of An et al. [32]. These results show that the BEM results are in a very good agreement to the FEM and NMM results. Thus, the validity of the proposed technique was confirmed.

6. Conclusions

Based on an explicit staggered predictor-corrector scheme The primary goal of this paper is to present a novel fractional-order theory that will aid in the advancement of functionally graded magnetic thermoelectric materials’ technological and industrial applications. Three-temperature nonlinear generalized thermoelasticity of functionally graded magnetic thermoelectric materials is the name given to this theory (FGMTMs). We must successfully adopt computerized numerical methods for solving and simulating complex nonlinear FGM problems in order to successfully guide the current research field toward the development of new functionally graded materials (FGMs). The governing equations are extremely difficult to solve experimentally or analytically due to the proposed theory’s severe nonlinearity. To address this issue, new numerical approaches for solving such equations must be developed. We propose a new formulation of the boundary element method for solving the theory’s governing equations (BEM). Because of the advantages of the BEM approach, such as the ability to deal with issues with complicated shapes that are difficult to deal with using standard methods, and the lack of the need for internal domain discretization. It also necessitates low CPU utilization and memory storage. As a result, the BEM is appropriate for a wide variety of sophisticated FGMTMs applications. The numerical results are discussed, with an emphasis on the effects of magnetic fields and graded parameters on thermal stresses in FGMTMs. Based on the findings, the suggested BEM technique is more effective, precise, and stable than GFDM or FEM–NMM. Current numerical results for our problem may be of interest to computer scientists, material science researchers, engineers, and designers and developers of functionally graded materials..