A New BEM for Modeling and Optimization of 3T Fractional Nonlinear Generalized Magneto-Thermoelastic Multi-Material ISMFGA Structures Subjected to Moving Heat Source

The main purpose of this chapter, which represents one of the chapters of a fractal analysis book, is to propose a new boundary element method (BEM) formulation based on time fractional order theory of thermoelasticity for modeling and optimization of three temperature (3T) multi-material initially stressed multilayered functionally graded anisotropic (ISMFGA) structures subjected to moving heat source. Fractional order derivative considered in the current chapter has been found to be an accurate mathematical tool for solving the difficulty of our physical and numerical modeling. Furthermore, this chapter shed light on the practical application aspects of boundary element method analysis and topology optimization of fractional order thermoelastic ISMFGA structures. Numerical examples based on the multi-material topology optimization algorithm and bi-evolutionary structural optimization method (BESO) are presented to study the effects of fractional order parameter on the optimal design of thermoelastic ISMFGA structures. The numerical results are depicted graphically to show the effects of fractional order parameter on the sensitivities of total temperature, displacement components and thermal stress components. The numerical results also show the effects of fractional order parameter on the final topology of the ISMFGA structures and demonstrate the validity and accuracy of our proposed technique.


Introduction
The fractional calculus has recently been widely used to study the theory and applications of derivatives and integrals of arbitrary non-integer order. This branch of mathematical analysis has emerged in recent years as an effective and powerful tool for the mathematical modeling of various engineering, industrial, and materials science applications [1][2][3]. The fractional-order operators are useful in describing the memory and hereditary properties of various materials and processes, due to their nonlocal nature. It clearly reflects from the related literature produced by leading fractional calculus journals that the primary focus of the investigation had shifted from classical integer-order models to fractional order models [4,5]. Fractional calculus has important applications in hereditary solid mechanics, fluid dynamics, viscoelasticity, heat conduction modeling and identification, biology, food engineering, econophysics, biophysics, biochemistry, robotics and control theory, signal and image processing, electronics, electric circuits, wave propagation, nanotechnology, etc. [6][7][8].
Numerous mathematicians have contributed to the history of fractional calculus, where Euler mentioned interpolating between integral orders of a derivative in 1730. Then, Laplace defined a fractional derivative by means of an integral in 1812.
Lacroix introduced the first fractional order derivative which appeared in a calculus in 1819, where he expressed the nth derivative of the function y ¼ x m as follows: Liouville assumed that d v dx v e ax ð Þ ¼ a v e ax for v > 0 to obtain the following fractional order derivative: Laurent has been using the Cauchy's integral formula for complex valued analytical functions to define the integration of arbitrary order v > 0 as follows: where c D v x denotes differentiation of order v of the function f along the x-axis. Cauchy introduced the following fractional order derivative: Caputo introduced his fractional derivative of order α < 0 to be defined as follows: physics, nuclear reactors, and other industrial applications. Due to computational difficulties in solving nonlinear generalized magneto-thermoelastic problems in general analytically, many numerical techniques have been developed and implemented for solving such problems [9][10][11][12][13][14][15][16][17]. The boundary element method (BEM) [18][19][20][21][22][23][24][25][26][27][28][29][30][31] has been recognized as an attractive alternative numerical method to domain methods [32][33][34][35][36] like finite difference method (FDM), finite element method (FEM), and finite volume method (FVM) in engineering applications. The superior feature of BEM over domain methods is that only the boundary of the domain needs to be discretized, which often leads to fewer elements and easier to use. This advantage of BEM over other domain methods has significant importance for modeling and optimization of thermoelastic problems which can be implemented using BEM with little cost and less input data. Nowadays, the BEM has emerged as an accurate and efficient computational technique for solving complicated inhomogeneous and non-linear problems in physical and engineering applications .
In the present chapter, we introduce a practical engineering application of fractal analysis in the field of thermoelasticity, where the thermal field is described by time fractional three-temperature radiative heat conduction equations. Fractional order derivative considered in the current chapter has high ability to remove the difficulty of our numerical modeling. A new boundary element method for modeling and optimization of 3T fractional order nonlinear generalized thermoelastic multi-material initially stressed multilayered functionally graded anisotropic (ISMFGA) structures subjected to moving heat source is investigated. Numerical results show that the fractional order parameter has a significant effect on the sensitivities of displacements, total three-temperature, and thermal stresses. Numerical examples show that the fractional order parameter has a significant effect on the final topology of ISMFGA structures. Numerical results of the proposed model confirm the validity and accuracy of the proposed technique, and numerical examples results demonstrate the validity of the BESO multi-material topology optimization method.
A brief summary of the chapter is as follows: Section 1 introduces the background and provides the readers with the necessary information to books and articles for a better understanding of fractional order problems and their applications. Section 2 describes the physical modeling of fractional order problems in three-temperature nonlinear generalized magneto-thermoelastic ISMFGA structures. Section 3 outlines the BEM implementation for modeling of 3T fractional nonlinear generalized magneto-thermoelastic problems of multi-material ISMFGA structures subjected to moving heat source. Section 4 introduces an illustration of the mechanisms of solving design sensitivities and optimization problem of the current chapter. Section 5 presents the new numerical results that describe the effects of fractional order parameter on the problem's field variations and on the final topology of multi-material ISMFGA structures.

Formulation of the problem
Consider a multilayered structure with n functionally graded layers in the xy-plane of a Cartesian coordinate. The x-axis is the common normal to all layers as shown in Figure 1. The thickness of the layer is denoted by h. The considered multilayered structure has been placed in a primary magnetic field H 0 acting in the direction of the y-axis.
According to the three-temperature theory, the governing equations of nonlinear generalized magneto-thermoelasticity in an initially stressed multilayered functionally graded anisotropic (ISMFGA) structure for the ith layer can be written in the following form: According to Fahmy [10], the time fractional order two-dimensional threetemperature (2D-3 T) radiative heat conduction equations in nondimensionless form can be expressed as follows: where in which The total energy of unit mass can be described by where σ ab , τ ab , and u i k are mechanical stress tensor, Maxwell's electromagnetic stress tensor, and displacement vector in the ith layer, respectively, c α (α = c, I, p) are constant T i α0 , T i α , C i abfg , and β i ab are, respectively, reference temperature, temperature, constant elastic moduli, and stress-temperature coefficients in the ith layer: μ i ,h, P i , ρ i , and c i sα are, respectively, magnetic permeability, perturbed magnetic field, initial stress, density, isochore specific heat coefficients in the ith layer; τ is the time; τ 0 and τ 1 are the relaxation times; i ¼ 1, 2, … , n represents the parameters in multilayered structure; and m is a functionally graded parameter. Also, we considered in the current study that the medium is subjected to a moving heat source of constant strength moving along x-axis with a constant velocity v. This moving heat source is assumed to have the following form: where, Q 0 is the heat source strength and δ is the delta function.

BEM numerical implementation
By using Eqs. (7)-(9), we can write (6) as where inertia term, temperature gradient, and initial stress terms are treated as the body forces.
In this section, we are interested in using a boundary element method for modeling the two-dimensional three-temperature radiation heat conduction equations coupled with electron, ion, and phonon temperatures.
According to finite difference scheme of Caputo at times f þ 1 ð ÞΔt and f Δτ, we obtain [1].
Based on Eq. (17), the fractional order heat Eq. (10) can be replaced by the following system: where j ¼ 1, 2, … :, F, f ¼ 0, 1, 2, … , F. Now, according to Fahmy [10], and applying the fundamental solution which satisfies (19), the boundary integral equations corresponding to (10) without heat sources can be expressed as Thus, the governing equations can be written in operator form as follows: where the operators L gb , f gb , L ab , and f ab are as follows: The differential Eq. (21) can be solved using the weighted residual method (WRM) to obtain the following integral equation: Now, the fundamental solution u i * df and traction vectors t i * da and t i a can be written as follows: Using integration by parts and sifting property of the Dirac distribution for (26), then using Eqs. (27) and (29), we can write the following elastic integral representation formula: The fundamental solution T i * can be defined as By using WRM and integration by parts, we can write (23) as follows: where By the use of sifting property, we obtain from (32) the thermal integral representation formula By combining (30) and (35), we obtain The nonlinear generalized magneto-thermoelastic vectors can be written in contracted notation form as By using the above vectors, we can express (36) as The source vector S A can be divided as where The representation formula (36) can also be written in matrix form as follows: In order to convert the domain integral in (42) into the boundary, we approximate the source vector S A by a series of known functions f q AE and unknown coefficients α q E as Thus, the representation formula (42) can be written as follows: By implementing the WRM to the following equations Then the elastic and thermal representation formulae are given as follows (Fahmy [46]): The representation formulae (55) and (56) can be combined into the following single equation: By substituting from Eq. (57) into Eq. (52), we obtain the following BEM coupled thermoelasticity formula: In order to compute the displacement sensitivity, Eq. (58) is differentiated with respect to ξ l as follows: According to the procedure of Fahmy [44], we can write (58) in the following form: The generalized displacements and velocities are approximated in terms of known tensor functions f q FD and unknown coefficients γ q D andγ q D : where Now, the gradients of the generalized displacements and velocities can also be approximated in terms of the tensor function derivatives as By substituting (63) into Eq. (45), we get By applying the point collocation procedure of Gaul et al. [43] to Eqs. (51) and (61), we obtain Similarly, applying the same point collocation procedure to Eqs. (64), (46), (47), (48), and (49) yields S where ψ, Γ AF , δ AF , and Ⅎ are assembled using the submatrices ψ ½ , Γ AF ½ , δ AF ½ , and Ⅎ ½ , respectively. Solving the system (65) for α and γ yields Now, the coefficient α can be written in terms of the unknown displacements U i , velocities _ U i , and accelerations € U i as An implicit-implicit staggered algorithm has been implemented for use with the BEM to solve the governing equations which can now be written in a suitable form after substitution of Eq. (72) into Eq. (60) as In many applications, the coupling term  z}|{ € U i nþ1 that appear in the heat conduction equation is negligible. Therefore, it is easier to predict the temperature than the displacement.
Hence Eqs. (73) and (74) lead to the following coupled system of differentialalgebraic equations (DAEs): Þ is the predicted temperature. By integrating Eq. (73) and using Eq. (75), we get From Eq. (77) we have where Substituting Eq. (79) into Eq. (78), we derive Substituting _ U i nþ1 from Eq. (79) into Eq. (75), we obtain By integrating the heat Eq. (74) and using Eq. (76), we obtain From Eq. (82) we get Substituting Eq. (84) into Eq. (83), we obtain Substituting _ T i nþ1 from Eq. (84) into Eq. (76), we get Now, a displacement predicted staggered procedure for the solution of (80) and (85) is as follows: The first step is to predict the propagation of the displacement wave field:  (86), respectively. The continuity conditions for temperature, heat flux, displacement, and traction that have been considered in the current chapter can be expressed as where n is the total number of layers, t a are the tractions which is defined by t a ¼ σ ab n b , and i ¼ 1, 2, … , n À 1.
The initial and boundary conditions of the present study are

Design sensitivity and optimization
According to Fahmy [58,60], the design sensitivities of displacements components and total 3T can be performed by implicit differentiation of (75) and (76), respectively, which describe the structural response with respect to the design variables, and then we can compute thermal stresses sensitivities.
The bi-directional evolutionary structural optimization (BESO) is the evolutionary topology optimization method that allows modification of the structure by either adding or removing material to or from the structure design. This addition or removal depends on the sensitivity analysis. Sensitivity analysis is the estimation of the response of the structure to the modification of design variables and is dependent on the calculation of derivatives [70][71][72][73][74][75][76][77][78][79][80].
The homogenized vector of thermal expansion coefficients α H can be written in terms of the homogenized elasticity matrix D H and the homogenized vector of stress-temperature coefficients β H as follows: For the material design, the derivative of the homogenized vector of thermal expansion coefficients can be written as where and where Ω j j is the volume of the base cell.

Numerical examples, results, and discussion
In order to show the numerical results of this study, we consider a monoclinic graphite-epoxy as an anisotropic thermoelastic material which has the following physical constants [57].
The elasticity tensor is expressed as Mass density ρ ¼ 7820kg=m 3 and heat capacity c = 461 J/kg K. The proposed technique that has been utilized in the present chapter can be applicable to a wide range of three-temperature nonlinear generalized thermoelastic problems of ISMFGA structures. The main aim of this chapter was to assess the impact of fractional order parameter on the sensitivities of total three-temperature, displacement components, and thermal stress components. Figure 2 shows the variation of the total temperature sensitivity along the x-axis. It was shown from this figure that the fraction order parameter has great effects on the total three-temperature sensitivity. Figures 3 and 4 show the variation of the displacement components u 1 and u 2 along the x-axis for different values of fractional order parameter. It was noticed from these figures that the fractional order parameter has great effects on the displacement sensitivities. Figures 5-7 show the variation of the thermal stress components σ 11 , σ 12 , and σ 22 , respectively, along the x-axis for different values of fractional order parameter. It was noted from these figures that the fractional order parameter has great influences on the thermal stress sensitivities.
Since there are no available results for the three-temperature thermoelastic problems, except for Fahmy's research [10][11][12][13][14]. For comparison purposes with the special cases of other methods treated by other authors, we only considered one-dimensional numerical results of the considered problem. In the special case under consideration, Figure 2. Variation of the total 3T sensitivity along x-axis. the displacement u 1 and thermal stress σ 11 results are plotted in Figures 8 and 9. The validity and accuracy of our proposed BEM technique were demonstrated by comparing our BEM results with the FEM results of Xiong and Tian [81], it can be noticed that the BEM results are found to agree very well with the FEM results. Example 1. Short cantilever beam.
The mean compliance has been minimized, to obtain the maximum stiffness, when the structure is subjected to moving heat source. In this example, we consider a short cantilever beam shown in Figure 10, where the BESO final topology of considered short cantilever beam shown in Figure 11a for α ¼ 0:5 and shown in Figure 11b for α ¼ 1:0. It is noticed from this figure that the fractional order parameter has a significant effect on the final topology of the multi-material ISMFGA structure.

Example 2. MBB beam.
It is known that extraordinary thermo-mechanical properties can be accomplished by combining more than two materials phases with conventional materials [75]. For this reason, it is essential that the topology optimization strategy permits more than two materials phases within the design domain. In this example, we consider a MBB beam shown in Figure 12, where the BESO final topology of MBB beam has been shown in Figure 13a for α ¼ 0:5 and shown in Figure 13b for α ¼ 1:0 to show the effect of fractional order parameter on the final topology of the multimaterial ISMFGA structure. Example 3. Roller-supported beam.  In this example, we consider a roller-supported beam shown in Figure 14, where the BESO final topology of a roller-supported beam shown in Figure 15a for α ¼ 0:5 and shown in Figure 15b for α ¼ 1:0. Example 4. Cantilever beam (validation example). In order to demonstrate the validity of our implemented BESO topology optimization technique, we consider isotropic case of a cantilever beam shown in Figure 16 as a special case of our anisotropic study to interpolate the elasticity matrix and the stress-temperature coefficients using the design variables X M , then we compare our      BESO final topology shown in Figure 17a with the material interpolation scheme of the solid isotropic material with penalization (SIMP) shown in Figure 17b.
The BESO topology optimization problem implemented in Examples 1 and 4, to find the distribution of the M material phases, with the volume constraint can be stated as where X M is the design variable; C M is the mean compliance; P is the total load on the structure, which is the sum of mechanical and thermal loads; u M is the displacement vector; V M * , is the volume of the solid material; N is the total number of elements; K M is the global stiffness matrix; x min is a small value (e.g., 0.0001), which it guarantee that none of the elements will be removed completely from design domain; f M,mec is the mechanical load vector; and f M,ter is the thermal load vector. Also, the BESO parameters considered in Examples 1-4 can be seen in Tables 1-4, respectively.
The BESO topology optimization problem implemented in Examples 2 and 3, to find the distribution of the two materials in the design domain, which minimize the compliance of the structure, subject to a volume constraint in both phases can be stated as.
Find  where V M * , j is the volume of jth material phase and i and j denote the element ith which is made of jth material.

Conclusion
The main purpose of this chapter is to describe a new boundary element formulation for modeling and optimization of 3T time fractional order nonlinear generalized thermoelastic multi-material ISMFGA structures subjected to moving heat source, where we used the three-temperature nonlinear radiative heat conduction equations combined with electron, ion, and phonon temperatures.
Numerical results show the influence of fractional order parameter on the sensitivities of the study's fields. The validity of the present method is examined and demonstrated by comparing the obtained outcomes with those known in the literature. Because there are no available data to confirm the validity and accuracy of our proposed technique, we replace the three-temperature radiative heat conduction with one-temperature heat conduction as a special case from our current general study of three-temperature nonlinear generalized thermoelasticity. In the considered special case of 3T time fractional order nonlinear generalized thermoelastic multi-material ISMFGA structures, the BEM results have been compared graphically with the FEM results; it can be noticed that the BEM results are in excellent agreement with the FEM results. These results thus demonstrate the validity and accuracy of our proposed technique. Numerical examples are solved using the multi-material topology optimization algorithm based on the bi-evolutionary structural optimization method (BESO). Numerical results of these examples show that the fractional order parameter affects the final result of optimization. The implemented optimization algorithm has proven to be an appropriate computational tool for material design.
Nowadays, the knowledge of 3T fractional order optimization of multi-material ISMFGA structures, can be utilized by mechanical engineers for designing heat exchangers, semiconductor nano materials, thermoelastic actuators, shape memory actuators, bimetallic valves and boiler tubes. As well as for chemists to observe the chemical processes such as bond breaking and bond forming.