Fractional BEM Analysis of Nonlinear Three-Temperature Functionally Graded Magnetic Thermoelectric Materials

1. Introduction

Thermal transport in solids is closely linked to lattice vibrations. Initial descriptions of quantum mechanics that were developed for electromagnetic radiation and heat capacity of solids formed the framework for the study of quantized vibrational energies of crystals. These lattice vibrations, which are typically modeled by phonons, transport thermal and acoustic energies in solids and thus constitute an essential aspect of their thermal and mechanical response. This is particularly significant in multi-field materials, where the performance is contingent upon thermal, electrical, magnetic, and elastic interactions.

Many materials exhibit pronounced changes in their effective properties when subjected to electric or magnetic fields. This behavior has motivated the development of electrorheological and magnetically responsive fluids, electro-active elastomers, and magneto-active elastomers for applications such as actuators, clutches, shock absorbers, valves, and adaptive devices. In solids, the coupling between field-induced changes in microstructure and the resulting elastic response provides a useful mechanism for designing multifunctional materials with controllable stiffness and compliance.

Fractional calculus (FC) extends differentiation and integration to non-integer orders and has become a useful modeling framework for systems with memory and hereditary effects [4]. Fractional formulations have been applied to circuit design [5], vibration analysis [7], hydro-turbine dynamics [8], control engineering [9,10], nanotechnology [11], and biological processes [12,13]. Fractional adaptive algorithms have also been developed for parameter estimation [1,2,3,14,15,16], where the fractional order can influence convergence and estimation behavior. These developments motivate the use of fractional operators in coupled thermomechanical models whose transient response cannot be represented adequately by classical integer-order descriptions.

The foundations of magnetoelasticity were established by Knopoff [17], Chadwick [18], and Kaliski and Petykiewicz [19]. With the growth of high-performance computing, numerical simulation has become an increasingly important complement to magnetic experiments, particularly when direct testing is difficult, hazardous, or expensive. Accurate computational models are therefore essential for resolving coupled magnetic, thermal, and mechanical fields in advanced materials.

The interaction of magnetic, thermal, and deformation fields in thermoelastic solids is relevant to geophysics, plasma physics, and multifunctional material systems. In many formulations, the magnetic field enters the mechanical equations through Lorentz forces and the electrical response through a modified Ohm law. Generalized thermoelastic models extend classical heat conduction by incorporating relaxation and non-Fourier effects. Ezzat and Awad developed a micropolar generalized magneto-thermoelastic model based on modified Ohm and Fourier laws. Three-temperature theories have subsequently been investigated for micropolar thermoelasticity [21], carbon-nanotube-reinforced composites [22], micropolar piezothermoelasticity [23], micropolar magneto-thermoviscoelasticity [24], and magneto-thermoviscoelasticity [25]. The strong nonlinearity and fractional character of these coupled models generally preclude closed-form solutions, making robust numerical techniques essential.

The present study constitutes a comprehensive re-examination, reconstruction, and substantive extension of earlier work undertaken within the same research program. The theoretical formulation, analytical assumptions, computational procedures, interpretation of results, and overall presentation have been systematically reassessed and, where necessary, reformulated to establish a more general, internally consistent, and methodologically rigorous framework. Particular attention has been devoted to identifying and resolving technical inaccuracies, methodological inconsistencies, and interpretive limitations present in the earlier work. All revisions, corrections, and substantive modifications presented in the current article were undertaken solely by the author of this work. Importantly, the formulation developed herein is broader in scope than the earlier collaborative analysis, which may be recovered as a restricted special case of the present framework when the corresponding simplifying assumptions, parameter restrictions, and constitutive conditions are imposed. The current contribution therefore extends beyond a conventional correction or incremental refinement and instead provides a comprehensive theoretical and computational framework within which the previous formulation is naturally embedded as a limiting case. Its principal objective is to correct and strengthen the scientific record, improve methodological reliability and transparency, and provide a more general basis for the analysis of the class of problems considered. In this respect, the present work reflects the author's commitment to the accuracy, integrity, completeness, and continuing improvement of the scientific record.

This work develops a fractional-order formulation for functionally graded magnetic thermoelectric materials within a nonlinear three-temperature generalized thermoelastic framework. The principal numerical contribution is a boundary element treatment of the coupled governing equations. By reducing the spatial discretization burden on the domain and retaining the geometric flexibility of BEM, the formulation is well suited to complex functionally graded structures. The numerical study focuses on the influence of the fractional-order parameter and the functional grading parameter on the thermal-stress field and evaluates the proposed implementation through special-case comparisons with previously reported GFDM and FEM–NMM results.

2. Formulation of the Problem

Consider the two-dimensional functionally graded magnetic thermoelectric structure shown in Figure 1, with thickness $h$ and length $L$ in a Cartesian coordinate system. The structure is subjected to an electric potential $\Phi (x,z,t)$ and placed in a constant external magnetic field. The material occupies a region $\Omega$ bounded by $\Gamma$, whose boundary subsets are defined according to the prescribed thermal and mechanical conditions.

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)=\left\{\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}\right.$$

(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}} \left(7\right)$

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)

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

(9b)

$${\left.{\mathbb{K}}_{\alpha }{\displaystyle \frac{\partial T_{\alpha }}{\partial n}}\right|}_{{\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{\left.T_{\alpha }^i{T}_{\alpha }^{*}\right|}_{\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)

Upon using the integration by parts once more 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}=\left\{\begin{array}{c}{\widehat{\mathbb{H}}}^{ij} if i\ne j\\ {\widehat{\mathbb{H}}}^{ij}+C^i if i=j\end{array}\right.$$

(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)

The coupled systems of algebraic equations are coded into MATLAB (R2018a) using an explicit staggered predictor-corrector scheme. During each time step, the obtained linear systems of equations are solved by the communication-avoiding generalized minimal residual (CA-GMRES) approach in order to determine the temperature, displacement, and microrotation fields.

5. Numerical Results and Discussion

The numerical study examines the sensitivity of the three stress components in response to the variations of the fractional-order and functional grading parameters $\mathsf{\alpha }$ and $\mathrm{m}$. Furthermore, the performance features and the correctness of the model under particular conditions are studied. The BEM discretization is carried out by employing 42 boundary elements and 68 internal nodes, which are depicted in Figure 2. The following analysis will be made based on the trends observed on the considered x coordinate.

The effects of the fractional-order parameter for $\mathsf{\alpha }=\mathrm{0.3,0.6,1.0}$ are depicted in Figure 3, Figure 4 and Figure 5. As seen from Figure 3, ${\mathsf{\sigma }}_{11}$ stays slightly compressive in most part of the domain and then starts decreasing very fast near the $\mathrm{x}$ coordinate value approximately equal to 1.9 until recovering to the $\mathrm{x}$ coordinate value of 2. Increase of $\mathsf{\alpha }$ leads to more negative values of ${\mathsf{\sigma }}_{11}$, especially in the high-gradient zone.

For ${\sigma }_{12}$ (Figure 4), the stress decreases with increasing $x$ and approaches zero near $x=2$. At a fixed position, a larger $\alpha$ produces a larger positive ${\sigma }_{12}$ over most of the plotted interval. The separation between the curves is most pronounced away from the right boundary and narrows as the stress approaches its terminal value.

The ${\sigma }_{22}$ response in Figure 5 is compressive and increases toward zero as $x$ increases. In contrast to ${\sigma }_{12}$, increasing $\alpha$ moves ${\sigma }_{22}$ downward, i.e., toward larger compressive magnitude. The fractional-order parameter therefore affects the individual stress components differently rather than producing a uniform amplification or attenuation of the entire stress tensor.

Figure 6, Figure 7 and Figure 8 study the effect of the functional grading coefficient $\mathrm{m}=\mathrm{0.3,0.6,1.0}$. As can be seen in Figure 6, ${\mathsf{\sigma }}_{11}$ rises with an increase in $\mathrm{x}$, achieves its local maximum for $\mathrm{x}\approx 1.9$, and falls off toward the right edge of the figure. An increase in $\mathrm{m}$ not only enhances ${\mathsf{\sigma }}_{11}$ but also amplifies the value of the local maximum.

As can be seen from Figure 7, ${\sigma }_{12}$ is a negative quantity over the plotted interval and tends to zero when $x$ → 2. Increasing $m$ makes ${\sigma }_{12}$ more negative over most of the interval, while the three curves converge near the right boundary. The grading effect is therefore strongest in the interior and weakens close to the terminal boundary.

A similar grading sensitivity is observed for ${\sigma }_{22}$ in Figure 8. Stress is still compressive but gets less and less negative as $x$ is increasing; on the other hand, the bigger $m$ is, the more negative ${\sigma }_{22}$ value will be for a particular point. Contrary to ${\sigma }_{12}$ , there is no confusion in distinguishing between the curves of  ${\sigma }_{22}$ .

When put together, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 illustrate the dependence of the thermoelastic deformation on the material property involved. An increase in $\alpha$ reduces ${\sigma }_{11}$ and ${\sigma }_{22}$ and increases ${\sigma }_{12}$, whereas an increase in $m$ increases ${\sigma }_{11}$ and reduces ${\sigma }_{12}$ and ${\sigma }_{22}$. This shows that the memory effect and spatial variation of the material can serve different means of rearranging the thermal stress components.

Table 1 summarizes the reported computational-resource comparison for the present BEM implementation and the GFDM and FEM–NMM data used in the manuscript. The BEM calculation uses 68 nodes and 42 boundary elements, compared with 50,000 nodes and 20,000 elements for GFDM and 48,000 nodes and 18,000 elements for FEM–NMM. The reported CPU-time, memory, and disk-space metrics are also markedly smaller for BEM. These values are consistent with the principal computational advantage of a boundary-focused discretization. Because the table does not specify hardware, physical units, or normalization for the resource and accuracy metrics, the numerical values should be interpreted as implementation-specific reported metrics rather than hardware-independent performance benchmarks.

No independent published solution is available in the manuscript for the complete coupled fractional three-temperature problem; therefore, numerical consistency is examined through limiting cases that can be compared with the reported GFDM and FEM–NMM data. Figure 9 compares ${\sigma }_{11}$ for $\alpha =0.6$ and $m=0$, corresponding to a homogeneous case, while Figure 10 considers $\alpha =0$ and $m=0.6$. In both comparisons, the BEM curve closely follows the reference curves over the full $x$-range, including the steep stress-gradient region near $x\approx 1.9$. The visible deviations remain small relative to the overall stress variation. These special-case comparisons support the numerical consistency of the proposed BEM formulation, while the full coupled model would still benefit from additional independent benchmark or experimental validation.

6. Conclusions

This work presented a fractional boundary element formulation for nonlinear three temperature analysis of functionally graded magnetic thermoelectric materials. The proposed model integrates coupled thermal, electromagnetic, and mechanical effects within a unified computational framework and provides a practical approach for treating complex material behavior without requiring complete discretization of the interior domain.

The numerical results demonstrate that both the fractional order and the material gradation strongly affect the distribution and intensity of thermal stresses. These parameters therefore provide important mechanisms for controlling the thermomechanical response of functionally graded thermoelectric structures. The comparison with established numerical methods for representative special cases shows close agreement and supports the accuracy and reliability of the proposed boundary element procedure.

The developed formulation also offers computational advantages in terms of reduced discretization requirements, storage demand, and computational effort. These features make the method suitable for the simulation and design of advanced thermoelectric structures subjected to coupled thermal, magnetic, electric, and mechanical loading. The present framework can also provide a foundation for future extensions involving more complex material architectures, transient loading conditions, and multiphysics applications.