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

Abstract

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.

1. 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 $\mathrm{nth}$ derivative of the function $\mathrm{y}={\mathrm{x}}^{\mathrm{m}}$ as follows:

$$\frac{{\mathrm{d}}^{\mathrm{n}}}{{\mathrm{dx}}^{\mathrm{n}}}=\frac{\mathrm{\Gamma }\left(\mathrm{m}+1\right)}{\mathrm{\Gamma }\left(\mathrm{m}‐\mathrm{n}+1\right)}{\mathrm{x}}^{\mathrm{m}‐\mathrm{n}}$$E1

Liouville assumed that $\frac{{\mathrm{d}}^{\mathrm{v}}}{{\mathrm{dx}}^{\mathrm{v}}}\left({\mathrm{e}}^{\mathrm{ax}}\right)={\mathrm{a}}^{\mathrm{v}}{\mathrm{e}}^{\mathrm{ax}}\text{for}\mathrm{v}>0$ to obtain the following fractional order derivative:

$$\frac{{\mathrm{d}}^{\mathrm{v}}{\mathrm{x}}^{‐\mathrm{a}}}{{\mathrm{dx}}^{\mathrm{v}}}={\left(‐1\right)}^{\mathrm{v}}\frac{\mathrm{\Gamma }\left(\mathrm{a}+\mathrm{v}\right)}{\mathrm{\Gamma }\left(\mathrm{a}\right)}{\mathrm{x}}^{‐\mathrm{a}‐\mathrm{v}}$$E2

Laurent has been using the Cauchy’s integral formula for complex valued analytical functions to define the integration of arbitrary order $\mathrm{v}>0$ as follows:

$${}_{\mathrm{c}}{\mathrm{D}}_{\mathrm{x}}^{\mathrm{v}}\mathrm{f}\left(\mathrm{x}\right)={}_{\mathrm{c}}{\mathrm{D}}_{\mathrm{x}}^{\mathrm{m}‐\mathrm{\rho }}\mathrm{f}\left(\mathrm{x}\right)=\frac{{\mathrm{d}}^{\mathrm{m}}}{{\mathrm{dx}}^{\mathrm{m}}}\left[\frac{1}{\mathrm{\Gamma }\left(\mathrm{\rho }\right)}{\int }_{\mathrm{c}}^{\mathrm{x}}{\left(\mathrm{x}-\mathrm{t}\right)}^{\mathrm{\rho }-1}\mathrm{f}\left(\mathrm{t}\right)\mathrm{dt}\right],0<\mathrm{\rho }\le 1$$E3

where ${}_{\mathrm{c}}{\mathrm{D}}_{\mathrm{x}}^{\mathrm{v}}$ denotes differentiation of order $\mathrm{v}$ of the function $\mathrm{f}$ along the $\mathrm{x}‐\text{axis}$.

Cauchy introduced the following fractional order derivative:

$${\mathrm{f}}_{+}^{\left(\mathrm{\alpha }\right)}=\int \mathrm{f}\left(\mathrm{\tau }\right)\frac{{\left(\mathrm{t}‐\mathrm{\tau }\right)}^{\mathrm{\alpha }‐1}}{\mathrm{\Gamma }\left(‐\mathrm{\alpha }\right)}\mathrm{d\tau }$$E4

Caputo introduced his fractional derivative of order $\mathrm{\alpha }<0$ to be defined as follows:

$${\mathrm{D}}_{\ast }^{\mathrm{\alpha }}\mathrm{f}\left(\mathrm{t}\right)=\frac{1}{\mathrm{\Gamma }\left(\mathrm{m}‐\mathrm{\alpha }\right)}{\int }_0^{\mathrm{t}}\frac{{\mathrm{f}}^{\left(\mathrm{m}\right)}\left(\mathrm{\tau }\right)}{{\left(\mathrm{t}‐\mathrm{\tau }\right)}^{\mathrm{\alpha }+1‐\mathrm{m}}}\mathrm{d\tau },\mathrm{m}-1<\mathrm{\alpha }<\mathrm{m},\mathrm{\alpha }>0$$E5

Recently, research on nonlinear generalized magneto-thermoelastic problems has received wide attention due to its practical applications in various fields such as geomechanics, geophysics, petroleum and mineral prospecting, earthquake engineering, astronautics, oceanology, aeronautics, materials science, fiber-optic communication, fluid mechanics, automobile industries, aircraft, space vehicles, plasma 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 [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69].

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.

2. Formulation of the problem

Consider a multilayered structure with $n$ functionally graded layers in the $\mathit{xy}‐\text{plane}$ of a Cartesian coordinate. The $x‐\text{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‐\text{axis}$.

Figure 1.

Geometry of the considered problem.

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:

$${\sigma }_{\mathit{ab},b}+{\tau }_{\mathit{ab},b}-{\mathrm{\Gamma }}_{\mathit{ab}}={\rho }^i{\left(x+1\right)}^m{\ddot{u}}_a^i$$E6

$${\sigma }_{\mathit{ab}}={\left(x+1\right)}^m\left[C_{\mathit{\text{abfg}}}^i{u}_{f,g}^i-{\beta }_{\mathit{ab}}^i\left(T_{\alpha }^i-T_{\alpha 0}^i+{\tau }_1{\dot{T}}_{\alpha }^1\right)\right]$$E7

$${\tau }_{\mathit{ab}}={\mu }^i{\left(x+1\right)}^m\left({\tilde{h}}_a{H}_b+{\tilde{h}}_b{H}_a-{\delta }_{\mathit{ba}}\left({\tilde{h}}_f{H}_f\right)\right)$$E8

$${\mathrm{\Gamma }}_{\mathit{ab}}=P^i{\left(x+1\right)}^m\left(\frac{\partial u_a^i}{\partial x_b}-\frac{\partial u_b^i}{\partial x_a}\right)$$E9

According to Fahmy [10], the time fractional order two-dimensional three-temperature (2D-3 T) radiative heat conduction equations in nondimensionless form can be expressed as follows:

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

where

$$\overline{\mathbb{W}}\left(r,\tau \right)=\left\{\begin{array}{ll}-{\rho }^i{\mathbb{W}}_{\mathit{eI}}\left(T_e^i-T_I^i\right)-{\rho }^i{\mathbb{W}}_{\mathit{er}}\left(T_e^i-T_p^i\right)+\overline{\overline{\mathbb{W}}},& \alpha =e,{\delta }_1=1\\ {\rho }^i{\mathbb{W}}_{\mathit{eI}}\left(T_e^i-T_I^i\right)+\overline{\overline{\mathbb{W}}},& \alpha =I,{\delta }_1=1\\ {\rho }^i{\mathbb{W}}_{\mathit{er}}\left(T_e^i-T_p^i\right)+\overline{\overline{\mathbb{W}}},& \alpha =p,{\delta }_1=T_p^3\end{array}\right.$$E11

in which

$$\overline{\overline{\mathbb{W}}}\left(r,\tau \right)=-{\mathbb{K}}_{\alpha }^i{\dot{T}}_{\alpha ,\mathit{ab}}^i+{\beta }_{\mathit{ab}}^i{T}_{a0}^i{\tau }_0{\ddot{u}}_{a,b}^i+{\rho }^i{c}_{\mathit{s\alpha }}^i{\tau }_0{\ddot{T}}_{\alpha }^i-Q\left(x,\tau \right)$$E12

$${\mathbb{W}}_{\mathit{eI}}={\rho }^i{\mathbb{A}}_{\mathit{eI}}{T}_e^{i\left(-\frac{2}{3}\right)},{\mathbb{W}}_{\mathit{er}}={\rho }^i{\mathbb{A}}_{\mathit{er}}{T}_e^{i\left(-\frac{1}{2}\right)},{\mathbb{K}}_{\alpha }={\mathbb{A}}_{\alpha }{T}_{\alpha }^{i\left(\frac{5}{2}\right)},\alpha =e,I,{\mathbb{K}}_p={\mathbb{A}}_p{T}_p^{i\left(3+\mathbb{B}\right)}$$E13

The total energy of unit mass can be described by

$$P=P_e+P_I+p_p,P_e=c_{\mathit{e}}{T}_e^i,P_I=c_{\mathit{I}}{T}_I^i,P_p=\frac{1}{4}{c}_{\mathit{p}}{T}_p^{4i}$$E14

where ${\sigma }_{\mathit{ab}}$, ${\tau }_{\mathit{ab}}$, and $u_k^i$ are mechanical stress tensor, Maxwell’s electromagnetic stress tensor, and displacement vector in the ith layer, respectively, c_α(α = c, I, p) are constant $T_{\alpha 0}^i$, $T_{\alpha }^i$, $C_{\mathit{\text{abfg}}}^i$, and ${\beta }_{\mathit{ab}}^i$ are, respectively, reference temperature, temperature, constant elastic moduli, and stress-temperature coefficients in the ith layer: ${\mu }^i$, $\tilde{h}$, $P^i$, ${\rho }^i$, and $c_{\mathit{s\alpha }}^i$ are, respectively, magnetic permeability, perturbed magnetic field, initial stress, density, isochore specific heat coefficients in the ith layer; $\tau$ is the time; ${\tau }_0$ and ${\tau }_1$ are the relaxation times; $i=1,2,\dots ,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‐\text{axis}$ with a constant velocity $v$. This moving heat source is assumed to have the following form:

$$Q\left(x,\tau \right)=Q_0\delta \left(x-\mathit{v\tau }\right)$$E15

where, $Q_0$ is the heat source strength and $\delta$ is the delta function.

3. BEM numerical implementation

By using Eqs. (7)–(9), we can write (6) as

$$L_{\mathit{gb}}{u}_f^i={\rho }^i{\ddot{u}}_a^i-\left(D_a{T}_{\alpha }^i-P^i\left(\frac{\partial u_b^i}{\partial x_a}-\frac{\partial u_a^i}{\partial x_b}\right)\right)=f_{\mathit{gb}}$$E16

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 $\left(f+1\right)\mathrm{\Delta }t$ and $f\mathrm{\Delta }\tau$, we obtain [1].

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

where

$$W_{a,0}=\frac{{\left(\mathrm{\Delta }\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)$$E18

Based on Eq. (17), the fractional order heat Eq. (10) can be replaced by the following system:

$$\begin{array}{l}{W}_{a,0}{T}_{\alpha }^{i\left(f+1\right)}\left(r\right)-{\mathbb{K}}_{\alpha }\left(x\right)T_{\alpha ,\mathit{II}}^{i\left(f+1\right)}\left(r\right)-{\mathbb{K}}_{\alpha ,I}\left(x\right)T_{\alpha ,I}^{i\left(f+1\right)}\left(r\right)=W_{a,0}{T}_{\alpha }^{i\left(f\right)}\left(r\right)-{\mathbb{K}}_{\alpha }\left(x\right)T_{\alpha ,\mathit{II}}^{i\left(f\right)}\left(r\right)\\ -{\mathbb{K}}_{\alpha ,I}\left(x\right)T_{\alpha ,J}^{i\left(f\right)}\left(r\right)-\sum _{j=1}^f{W}_{a,j}\left(T_{\alpha }^{i\left(f+1-j\right)}\left(r\right)-T_{\alpha }^{i\left(f-j\right)}\left(r\right)\right)+{\overline{\mathbb{W}}}_m^{i\left(f+1\right)}\left(x,\tau \right)+{\overline{\mathbb{W}}}_m^{i\left(f\right)}\left(x,\tau \right)\end{array}$$E19

where $j=1,2,\dots .,F,f=0,1,2,\dots ,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

$$T_{\alpha }^i\left(\xi \right)={\int }_S\left[T_{\alpha }^i{q}^{i\ast }-T_{\alpha }^{i\ast }{q}^i\right]\mathit{dC}-{\int }_R{f}_{\mathit{ab}}{T}_{\alpha }^{i\ast }\mathit{dR}$$E20

Thus, the governing equations can be written in operator form as follows:

$$L_{\mathit{gb}}{u}_f^i=f_{\mathit{gb}},$$E21

$$L_{\mathit{ab}}{T}_{\alpha }^i=f_{\mathit{ab}}$$E22

where the operators $L_{\mathit{gb}}$, $f_{\mathit{gb}}$, $L_{\mathit{ab}}$, and $f_{\mathit{ab}}$ are as follows:

$$L_{\mathit{gb}}{D}_{\mathit{abf}}\frac{\partial }{\partial x_b}+D_{\mathit{af}}+\mathrm{\Lambda }{D}_{a1f},L_{\mathit{ab}}=D_{\tau }^a$$E23

$$f_{\mathit{gb}}={\rho }^i{\ddot{u}}_a^i-\left(D_a{T}_{\alpha }^i-P^i\left(\frac{\partial u_b^i}{\partial x_a}-\frac{\partial u_a^i}{\partial x_b}\right)\right)$$E24

$$f_{\mathit{ab}}=\frac{{\mathbb{K}}_{\alpha }}{D}{T}_{\alpha }^i\frac{\partial }{\partial \tau }$$E25

where

$$\begin{array}{c}{D}_{\mathit{abf}}=C_{\mathit{\text{abfg}}}\epsilon ,\epsilon =\frac{\partial }{\partial x_g},D_{\mathit{af}}=\mu H_0^2\left(\frac{\partial }{\partial x_a}+{\delta }_{a1}\mathrm{\Lambda }\right)\frac{\partial }{\partial x_f},\\ D_a=-{\beta }_{\mathit{ab}}^i\left(\frac{\partial }{\partial x_b}+{\delta }_{b1}\mathrm{\Lambda }+{\tau }_1\left(\frac{\partial }{\partial x_b}+\mathrm{\Lambda }\right)\frac{\partial }{\partial \tau }\right),\mathrm{\Lambda }=\frac{m}{x+1}.\end{array}$$

The differential Eq. (21) can be solved using the weighted residual method (WRM) to obtain the following integral equation:

$$\underset{R}{\int }\left(L_{\mathit{gb}}{u}_f^i-f_{\mathit{gb}}\right)u_{\mathit{da}}^{i\ast }\mathit{dR}=0$$E26

Now, the fundamental solution $u_{\mathit{df}}^{i\ast }$ and traction vectors $t_{\mathit{da}}^{i\ast }$ and $t_a^i$ can be written as follows:

$$L_{\mathit{gb}}{u}_{\mathit{df}}^{i\ast }=-{\delta }_{\mathit{ad}}\delta \left(x,\xi \right)$$E27

$$t_{\mathit{da}}^{i\ast }=C_{\mathit{\text{abfg}}}{u}_{\mathit{df},g}^{i\ast }{n}_b$$E28

$$t_a^i=\frac{{\overline{t}}_a^i}{{\left(x+1\right)}^m}=\left(C_{\mathit{\text{abfg}}}{u}_{f,g}^i-{\beta }_{\mathit{ab}}^i\left(T_{\alpha }^i+{\tau }_1{T}_{\alpha }^i\right)\right)n_b$$E29

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:

$$u_d^i\left(\xi \right)=\underset{C}{\int }\left(u_{\mathit{da}}^{i\ast }{t}_a^i-t_{\mathit{da}}^{i\ast }{u}_a^i+u_{\mathit{da}}^{i\ast }{\beta }_{\mathit{ab}}^i{T}_{\alpha }^i{n}_b\right)\mathit{dC}-\underset{R}{\int }{f}_{\mathit{gb}}{u}_{\mathit{da}}^{i\ast }\mathit{dR}$$E30

The fundamental solution $T^{i\ast }$ can be defined as

$$L_{\mathit{ab}}{T}_{\alpha }^{i\ast }=-\delta \left(x,\xi \right)$$E31

By using WRM and integration by parts, we can write (23) as follows:

$$\underset{R}{\int }\left(L_{\mathit{ab}}{T}_{\alpha }^i{T}_{\alpha }^{i\ast }-L_{\mathit{ab}}{T}_{\alpha }^{i\ast }{T}_{\alpha }^i\right)\mathit{dR}=\underset{C}{\int }\left(q^{i\ast }{T}_{\alpha }^i-q^i{T}_{\alpha }^{i\ast }\right)\mathit{dC}$$E32

where

$$q^i=-{\mathbb{K}}_{\alpha }^i{T}_{\alpha ,b}^i{n}_a$$E33

$$q^{i\ast }=-{\mathbb{K}}_{\alpha }^i{T}_{\alpha ,b}^{i\ast }{n}_a$$E34

By the use of sifting property, we obtain from (32) the thermal integral representation formula

$$T_{\alpha }^i\left(\xi \right)=\underset{C}{\int }\left(q^{i\ast }{T}_{\alpha }^i-q^i{T}_{\alpha }^{i\ast }\right)\mathit{dC}-\underset{C}{\int }{f}_{\mathit{ab}}{T}_{\alpha }^{i\ast }\mathit{dR}$$E35

By combining (30) and (35), we obtain

$$\left[\begin{array}{c}{u}_d^i\left(\xi \right)\\ T_{\alpha }^i\left(\xi \right)\end{array}\right]=\underset{C}{\int }\left\{-\left[\begin{array}{cc}{t}_{\mathit{da}}^{i\ast }& -u_{\mathit{da}}^{i\ast }{\beta }_{\mathit{ab}}{n}_b\\ 0& -q^{i\ast }\end{array}\right]\left[\begin{array}{c}{u}_a^i\\ T_{\alpha }^i\end{array}\right]+\left[\begin{array}{cc}{u}_{\mathit{da}}^{i\ast }& 0\\ 0& -T_{\alpha }^{i\ast }\end{array}\right]\left[\begin{array}{c}{t}_a^i\\ q^i\end{array}\right]\right\}\mathit{dC}-\underset{R}{\int }\left[\begin{array}{cc}{u}_{\mathit{da}}^{i\ast }& 0\\ 0& -T_{\alpha }^{i\ast }\end{array}\right]\left[\begin{array}{c}{f}_{\mathit{gb}}\\ -f_{\mathit{ab}}\end{array}\right]\mathit{dR}$$E36

The nonlinear generalized magneto-thermoelastic vectors can be written in contracted notation form as

$$U_A^i=\left\{\begin{array}{cc}{u}_a^i& a=A=1,2,3\\ T_{\alpha }^i& A=4\end{array}\right.$$E37

$${\mathrm{T}}_{\mathit{\alpha A}}^i=\left\{\begin{array}{cc}{t}_a^i& a=A=1,2,3\\ q^i& A=4\end{array}\right.$$E38

$$U_{\mathit{DA}}^{i\ast }=\left\{\begin{array}{ll}{u}_{\mathit{da}}^{i\ast }& d=D=1,2,3;a=A=1,2,3\\ 0& d=D=1,2,3;A=4\\ 0& D=4;a=A=1,2,3\\ -T_{\alpha }^{i\ast }& D=4;A=4\end{array}\right.$$E39

$${\tilde{T}}_{\mathit{\alpha DA}}^{i\ast }=\left\{\begin{array}{ll}{t}_{\mathit{da}}^{i\ast }& d=D=1,2,3;a=A=1,2,3\\ -{\tilde{u}}_d^{i\ast }& d=D=1,2,3;A=4\\ 0& D=4;a=A=1,2,3\\ -q^{i\ast }& D=4;A=4\end{array}\right.$$E40

$${\tilde{u}}_d^{i\ast }=u_{\mathit{da}}^{i\ast }{\beta }_{\mathit{af}}{n}_f$$E41

By using the above vectors, we can express (36) as

$$U_D^i\left(\xi \right)=\underset{C}{\int }\left(U_{\mathit{DA}}^{i\ast }{\mathrm{T}}_{\mathit{\alpha A}}^i-{\tilde{T}}_{\mathit{\alpha DA}}^i{U}_A^i\right)\mathit{dC}-\underset{R}{\int }{U}_{\mathit{DA}}^{i\ast }{S}_A\mathit{dR}$$E42

The source vector $S_A$ can be divided as

$$S_A=S_A^0+S_A^T+S_A^u+S_A^{\dot{T}}+S_A^{\ddot{T}}+S_A^{\ddot{u}}$$E43

where

$$S_A^0=\left\{\begin{array}{ll}0& A=1,2,3\\ Q_0\delta \left(x-\mathit{v\tau }\right)& A=4\end{array}\right.$$E44

$$S_A^T={\omega }_{\mathit{AF}}{U}_F^i\text{with}{\omega }_{\mathit{AF}}=\left\{\begin{array}{ll}-D_a& A=1,2,3;F=4\\ \xi \text{∇}\left[{\mathbb{K}}_{\alpha }^i\text{∇}\right]& \text{otherwise}\end{array}\right.$$E45

$$S_A^u=\psi U_F^i\text{with}\psi =\left\{\begin{array}{ll}{P}^i\left(\frac{\partial }{\partial x_b}-\frac{\partial }{\partial x_a}\right)& A=1,2,3;F=1,2,3,\\ 0& A=4;F=4\end{array}\right.$$E46

$$S_A^{\dot{T}}={\mathrm{\Gamma }}_{\mathit{AF}}{\dot{U}}_F^i\text{with}{\mathrm{\Gamma }}_{\mathit{AF}}=\left\{\begin{array}{ll}-{\beta }_{\mathit{ab}}^i{\tau }_1\left(\frac{\partial }{\partial x_b}+\mathrm{\Lambda }\right)\frac{\partial }{\partial \tau }& A=4;F=4\\ -{\mathbb{K}}_{\alpha }^i& \text{otherwise}\end{array}\right.$$E47

$$S_A^{\ddot{T}}={\delta }_{\mathit{AF}}{\ddot{U}}_F^i\text{with}{\delta }_{\mathit{AF}}=\left\{\begin{array}{ll}0& A=4;F=4\\ {\rho }^i{c}_{\mathit{s\alpha }}^i{\mathit{\tau }}_0& \text{otherwise}\end{array}\right.$$E48

$$S_A^{\ddot{u}}=\mathcal{Ⅎ}{\ddot{U}}_F^i\text{with}\mathcal{Ⅎ}=\left\{\begin{array}{ll}{\rho }^i& A=1,2,3;F=1,2,3,\\ {\beta }_{\mathit{ab}}^i{T}_{\alpha 0}^i{\tau }_0& A=4;F=4\end{array}\right.$$E49

The representation formula (36) can also be written in matrix form as follows:

$$\begin{array}{ll}\left[S_A\right]=& -\left[\begin{array}{c}0\\ Q_0\delta \left(x-\mathit{v\tau }\right)\end{array}\right]+\left[\begin{array}{c}-D_a{T}_{\alpha }^i\\ \xi \text{∇}\left[{\mathbb{K}}_{\alpha }^i\text{∇}{T}_{\alpha }^i\left(r,\tau \right)\right]\end{array}\right]+\left[\begin{array}{c}{P}^i\left(u_{b,a}^i-u_{a,b}^i\right)\\ 0\end{array}\right]\\ & +\left[\begin{array}{c}-{\beta }_{\mathit{ab}}^i{\tau }_1\left(\frac{\partial }{\partial x_b}+\mathrm{\Lambda }\right){\dot{T}}_{\alpha }^i\\ -{\mathbb{K}}_{\alpha }^i{\dot{T}}_{\alpha }^i\end{array}\right]+{\rho }^i{c}_{\mathit{s\alpha }}^i{\tau }_0\left[\begin{array}{c}0\\ {\ddot{T}}_{\alpha }^i\end{array}\right]+\left[\begin{array}{c}{\rho }^i{\ddot{u}}_a^i\\ {\beta }_{\mathit{ab}}^i{T}_{\alpha 0}^i{\tau }_0{\ddot{u}}_{f,g}^i\end{array}\right]\end{array}$$E50

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_{\mathit{AE}}^q$ and unknown coefficients ${\alpha }_E^q$ as

$$S_A\approx \sum _{q=1}^E{f}_{\mathit{AE}}^q{\alpha }_E^q$$E51

Thus, the representation formula (42) can be written as follows:

$$U_D\left(\xi \right)=\underset{C}{\int }\left(U_{\mathit{DA}}^{i\ast }{T}_{\mathit{\alpha A}}^i-{\tilde{T}}_{\mathit{\alpha DA}}^{i\ast }{U}_A^i\right)\mathit{dC}-\sum _{q=1}^N\underset{R}{\int }{U}_{\mathit{DA}}^{i\ast }{f}_{\mathit{AE}}^q\mathit{dR}{\alpha }_E^q$$E52

By implementing the WRM to the following equations

$$L_{\mathit{gb}}{u}_{\mathit{fe}}^{\mathit{iq}}=f_{\mathit{ae}}^q$$E53

$$L_{\mathit{ab}}{T}_{\alpha }^{\mathit{iq}}=f_{\mathit{pj}}^q$$E54

Then the elastic and thermal representation formulae are given as follows (Fahmy [46]):

$$u_{\mathit{de}}^{\mathit{iq}}\left(\xi \right)=\underset{C}{\int }\left(u_{\mathit{da}}^{i\ast }{t}_{\mathit{ae}}^{\mathit{iq}}-t_{\mathit{da}}^{i\ast }{u}_{\mathit{ae}}^{\mathit{iq}}\right)\mathit{dC}-\underset{R}{\int }{u}_{\mathit{da}}^{i\ast }{f}_{\mathit{ae}}^q\mathit{dR}$$E55

$$T_{\alpha }^{\mathit{iq}}\left(\xi \right)=\underset{C}{\int }\left(q^{i\ast }{T}_{\alpha }^{\mathit{iq}}-q^{\mathit{iq}}{T}_{\alpha }^{i\ast }\right)\mathit{dC}-\underset{R}{\int }{f}^q{T}_{\alpha }^{i\ast }\mathit{dR}$$E56

The representation formulae (55) and (56) can be combined into the following single equation:

$$U_{\mathit{DE}}^{\mathit{iq}}\left(\xi \right)=\underset{C}{\int }\left(U_{\mathit{DA}}^{i\ast }{T}_{\mathit{\alpha AE}}^{\mathit{iq}}-T_{\mathit{\alpha DA}}^{i\ast }{U}_{\mathit{AE}}^{\mathit{iq}}\right)\mathit{dC}-\underset{R}{\int }{U}_{\mathit{DA}}^{i\ast }{f}_{\mathit{AE}}^{\mathit{iq}}\mathit{dR}$$E57

By substituting from Eq. (57) into Eq. (52), we obtain the following BEM coupled thermoelasticity formula:

$$U_D^i\left(\xi \right)=\underset{C}{\int }\left(U_{\mathit{DA}}^{i\ast }{T}_{\mathit{\alpha A}}^i-{\stackrel{⌣}{T}}_{\mathit{\alpha DA}}^{i\ast }{U}_A^i\right)\mathit{dC}+\sum _{q=1}^E\left(U_{\mathit{DE}}^{\mathit{iq}}\left(\xi \right)+\underset{C}{\int }\left(T_{\mathit{\alpha DA}}^{i\ast }{U}_{\mathit{AE}}^{\mathit{iq}}-U_{\mathit{DA}}^{i\ast }{T}_{\mathit{\alpha AE}}^{\mathit{iq}}\right)\mathit{dC}\right){\alpha }_E^q$$E58

In order to compute the displacement sensitivity, Eq. (58) is differentiated with respect to ${\xi }_l$ as follows:

$$\frac{\partial U_D^i\left(\xi \right)}{\partial {\xi }_l}=-\underset{C}{\int }\left(U_{\mathit{DA},l}^{i\ast }{T}_{\mathit{\alpha A}}^i-{\stackrel{⌣}{T}}_{\mathit{\alpha DA},l}^{i\ast }{U}_A^i\right)\mathit{dC}+\sum _{q=1}^E\left(\frac{\partial U_{\mathit{DE}}^{\mathit{iq}}\left(\xi \right)}{\partial {\xi }_l}-\underset{C}{\int }\left(T_{\mathit{\alpha DA},l}^{i\ast }{U}_{\mathit{\alpha AE}}^{\mathit{iq}}-U_{\mathit{DA},l}^{i\ast }{T}_{\mathit{\alpha AE}}^{\mathit{iq}}\right)\mathit{dC}\right){\alpha }_E^q$$E59

According to the procedure of Fahmy [44], we can write (58) in the following form:

$$\check{\zeta }U-\eta T_{\alpha }=\left(\zeta \check{U}-\eta \check{\wp }\right)\overline{\alpha }$$E60

The generalized displacements and velocities are approximated in terms of known tensor functions $f_{\mathit{FD}}^q$ and unknown coefficients ${\gamma }_D^q$ and ${\tilde{\gamma }}_D^q$:

$$U_F^i\approx \sum _{q=1}^N{f}_{\mathit{FD}}^q\left(x\right){\gamma }_D^q$$E61

where

$$f_{\mathit{FD}}^q=\left\{\begin{array}{cc}{f}_{\mathit{fd}}^q& f=F=1,2,3;d=D=1,2,3\\ f^q& F=4;D=4\\ 0& \text{otherwise}\end{array}\right.$$E62

Now, the gradients of the generalized displacements and velocities can also be approximated in terms of the tensor function derivatives as

$$U_{F,g}^i\approx \sum _{q=1}^N{f}_{\mathit{FD},g}^q\left(x\right){\gamma }_K^q$$E63

By substituting (63) into Eq. (45), we get

$$S_A^T=\sum _{q=1}^N{S}_{\mathit{AF}}{f}_{\mathit{FD},g}^q{\gamma }_D^q$$E64

By applying the point collocation procedure of Gaul et al. [43] to Eqs. (51) and (61), we obtain

$$\check{S}=J\overline{\alpha },U^i=J^{\prime }\gamma ,$$E65

Similarly, applying the same point collocation procedure to Eqs. (64), (46), (47), (48), and (49) yields

$${\check{S}}^{T_{\alpha }^i}={\mathcal{B}}^T\gamma$$E66

$$S_A^u=\psi U^i$$E67

$${\check{S}}^{T_{\alpha }^{\dot{\iota }}}={\overline{\mathrm{\Gamma }}}_{\mathit{AF}}{\dot{U}}^i$$E68

$${\check{S}}^{T_{\alpha }^{\ddot{\iota }}}={\overline{\delta }}_{\mathit{AF}}{\ddot{U}}^i$$E69

$${\check{S}}^{\ddot{u}}=\overline{\mathcal{Ⅎ}}{\ddot{U}}^i$$E70

where $\overline{\psi }$, ${\overline{\mathrm{\Gamma }}}_{\mathit{AF}}$, ${\overline{\delta }}_{\mathit{AF}}$, and $\overline{\mathcal{Ⅎ}}$ are assembled using the submatrices $\left[\psi \right]$, $\left[{\mathrm{\Gamma }}_{\mathit{AF}}\right]$, $\left[{\delta }_{\mathit{AF}}\right]$, and $\left[\mathcal{Ⅎ}\right]$, respectively.

Solving the system (65) for $\overline{\alpha }$ and $\gamma$ yields

$$\overline{\alpha }=J^{-1}\check{S},\gamma ={J^{\prime }}^{-1}{U}^i$$E71

Now, the coefficient $\overline{\alpha }$ can be written in terms of the unknown displacements $U^i$, velocities ${\dot{U}}^i$, and accelerations ${\ddot{U}}^i$ as

$$\overline{\alpha }=J^{-1}\left({\check{S}}^0+\left({\mathcal{B}}^T{J^{\prime }}^{-1}+\overline{\psi }\right)U^i+{\overline{\mathrm{\Gamma }}}_{\mathit{AF}}{\dot{U}}^i+\left(\overline{\mathcal{Ⅎ}}+{\overline{\delta }}_{\mathit{AF}}\right){\ddot{U}}^i\right)$$E72

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

$$\overbrace{M}{\ddot{U}}^i+\overbrace{\mathrm{\Gamma }}{\dot{U}}^i+\overbrace{K}{U}^i={\overbrace{\mathbb{Q}}}^i$$E73

$$\overbrace{X}{\ddot{T}}_{\alpha }^i+\overbrace{A}{\dot{T}}_{\alpha }^i+\overbrace{B}{T}_{\alpha }^i=\overbrace{\mathbb{Z}}{\ddot{U}}^i+\overbrace{\mathbb{R}}$$E74

where $V=\left(\eta \check{\wp }-\zeta \check{U}\right)J^{-1},$$\overbrace{M}=V\left(\overline{\mathcal{Ⅎ}}+{\overline{\delta }}_{\mathit{AF}}\right),$$\overbrace{\mathrm{\Gamma }}=V{\overline{\mathrm{\Gamma }}}_{\mathit{AF}},$$\overbrace{K}=-\check{\zeta }+V\left({\mathcal{B}}^T{J^{\prime }}^{-1}+\overline{\psi }\right)$, ${\overbrace{\mathbb{Q}}}^i=-\eta \check{T}+V{\check{S}}^0,$$\overbrace{X}=-{\rho }^i{c}_{\mathit{s\alpha }}^i{\tau }_0,$$\overbrace{A}=-{\mathbb{K}}_{\alpha }^i,$$\overbrace{B}=\xi \text{∇}\left[{\mathbb{K}}_{\alpha }^i\text{∇}\right],$$\overbrace{\mathbb{Z}}={\beta }_{\mathit{ab}}^i{T}_{\alpha 0}^i{\tau }_0,$$\overbrace{\mathbb{R}}=-Q_0\delta \left(x-\mathit{v\tau }\right)$.

where ${\ddot{U}}^i,{\dot{U}}^i,U^i,T^i$ and ${\overbrace{\mathbb{Q}}}^i$ are, respectively, acceleration, velocity, displacement, temperature, and external force vectors, and $V$,$\overbrace{M}$, $\overbrace{\mathrm{\Gamma }}$, $\overbrace{K}$, $\overbrace{\mathrm{A}}$, and $\overbrace{\mathrm{B}}$ are, respectively, volume, mass, damping, stiffness, capacity, and conductivity matrices.

In many applications, the coupling term $\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}^i$ 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 differential-algebraic equations (DAEs):

$$\overbrace{M}{\ddot{U}}_{n+1}^i+\overbrace{\mathrm{\Gamma }}{\dot{U}}_{n+1}^i+\overbrace{K}{U}_{n+1}^i={\overbrace{\mathbb{Q}}}_{n+1}^{\mathit{ip}}$$E75

$$\overbrace{\mathrm{X}}{\ddot{T}}_{\alpha \left(n+1\right)}^i+\overbrace{\mathrm{A}}{\dot{T}}_{\alpha \left(n+1\right)}^i+\overbrace{\mathrm{B}}{T}_{\alpha \left(n+1\right)}^i=\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}^i+\overbrace{\mathbb{R}}$$E76

where ${\overbrace{\mathbb{Q}}}_{n+1}^{\mathit{ip}}=\eta T_{\alpha \left(n+1\right)}^{\mathit{ip}}+V{\check{S}}^0$ and $T_{\alpha \left(n+1\right)}^{\mathit{ip}}$ is the predicted temperature.

By integrating Eq. (73) and using Eq. (75), we get

$$\begin{array}{ll}{\dot{U}}_{n+1}^i& ={\dot{U}}_n^i+\frac{\mathrm{\Delta }\tau }{2}\left({\ddot{U}}_{n+1}^i+{\ddot{U}}_n^i\right)\\ & ={\dot{U}}_n^i+\frac{\mathrm{\Delta }\tau }{2}\left[{\ddot{U}}_n^i+{\overbrace{M}}^{-1}\left({\overbrace{\mathbb{Q}}}_{n+1}^{\mathit{ip}}-\overbrace{\mathrm{\Gamma }}{\dot{U}}_{n+1}^i-\overbrace{K}{U}_{n+1}^i\right)\right]\end{array}$$E77

$$\begin{array}{ll}{U}_{n+1}^i& =U_n^i+\frac{\mathrm{\Delta }\tau }{2}\left({\dot{U}}_{n+1}^i+{\dot{U}}_n^i\right)\\ & =U_n^i+\mathrm{\Delta }\tau {\dot{U}}_n^i+\frac{\mathrm{\Delta }{\tau }^2}{4}\left[{\ddot{U}}_n^i+{\overbrace{M}}^{-1}\left({\overbrace{\mathbb{Q}}}_{n+1}^{\mathit{ip}}-\overbrace{\mathrm{\Gamma }}{\dot{U}}_{n+1}^i-\overbrace{K}{U}_{n+1}^i\right)\right]\end{array}$$E78

From Eq. (77) we have

$${\dot{U}}_{n+1}^i={\overline{\mathrm{\gamma }}}^{-1}\left[{\dot{U}}_n^i+\frac{\mathrm{\Delta }\tau }{2}\left[{\ddot{U}}_n^i+{\overbrace{M}}^{-1}\left({\overbrace{\mathbb{Q}}}_{n+1}^{\mathit{ip}}-\overbrace{K}{U}_{n+1}^i\right)\right]\right]$$E79

where $\overline{\mathrm{\gamma }}=\left(I+\frac{\mathrm{\Delta }\tau }{2}{\overbrace{M}}^{-1}\overbrace{\mathrm{\Gamma }}\right)$.

Substituting Eq. (79) into Eq. (78), we derive

$$\begin{array}{l}{U}_{n+1}^i=U_n^i+\mathrm{\Delta }\tau {\dot{U}}_n^i+\frac{\mathrm{\Delta }{\tau }^2}{4}\\ \left[{\ddot{U}}_n^i+{\overbrace{M}}^{-1}\left({\overbrace{\mathbb{Q}}}_{n+1}^{\mathit{ip}}-\overbrace{\mathrm{\Gamma }}{\overline{\mathrm{\gamma }}}^{-1}\left[{\dot{U}}_n^i+\frac{\mathrm{\Delta }\tau }{2}\left[{\ddot{U}}_n^i+{\overbrace{M}}^{-1}\left({\overbrace{\mathbb{Q}}}_{n+1}^{\mathit{ip}}-\overbrace{K}{U}_{n+1}^i\right)\right]\right]-\overbrace{K}{U}_{n+1}^i\right)\right]\end{array}$$E80

Substituting ${\dot{U}}_{n+1}^i$ from Eq. (79) into Eq. (75), we obtain

$$\begin{array}{l}{\ddot{U}}_{n+1}^i={\overbrace{M}}^{-1}\\ \left[{\overbrace{\mathbb{Q}}}_{n+1}^{\mathit{ip}}-\overbrace{\mathrm{\Gamma }}\left[{\overline{\mathrm{\gamma }}}^{-1}\left[{\dot{U}}_n^i+\frac{\mathrm{\Delta }\tau }{2}\left[{\ddot{U}}_n^i+{\overbrace{M}}^{-1}\left({\overbrace{\mathbb{Q}}}_{n+1}^{\mathit{ip}}-\overbrace{K}{U}_{n+1}^i\right)\right]\right]\right]-\overbrace{K}{U}_{n+1}^i\right]\end{array}$$E81

By integrating the heat Eq. (74) and using Eq. (76), we obtain

$$\begin{array}{ll}{\dot{T}}_{\alpha \left(n+1\right)}^i& ={\dot{T}}_n^i+\frac{\mathrm{\Delta }\tau }{2}\left({\ddot{T}}_{\alpha \left(n+1\right)}^i+{\ddot{T}}_{\mathit{\alpha n}}^i\right)\\ & ={\dot{T}}_{\mathit{\alpha n}}^i+\frac{\mathrm{\Delta }\tau }{2}\left({\overbrace{\mathrm{X}}}^{-1}\left[\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}^i+\overbrace{\mathbb{R}}-\overbrace{\mathrm{A}}{\dot{T}}_{\alpha \left(n+1\right)}^i-\overbrace{\mathrm{B}}{T}_{\alpha \left(n+1\right)}^i\right]+{\ddot{T}}_{\mathit{\alpha n}}^i\right)\end{array}$$E82

$$\begin{array}{ll}{T}_{\alpha \left(n+1\right)}^i& =T_{\mathit{\alpha n}}^i+\frac{\mathrm{\Delta }\tau }{2}\left({\dot{T}}_{\alpha \left(n+1\right)}^i+{\dot{T}}_{\mathit{\alpha n}}^i\right)\\ & =T_{\mathit{\alpha n}}^i+\mathrm{\Delta }\tau {\dot{T}}_{\mathit{\alpha n}}^i+\frac{\mathrm{\Delta }{\tau }^2}{4}\left({\ddot{T}}_{\mathit{\alpha n}}^i+{\overbrace{\mathrm{X}}}^{-1}\left[\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}^i+\overbrace{\mathbb{R}}-\overbrace{\mathrm{A}}{\dot{T}}_{\alpha \left(n+1\right)}^i-\overbrace{\mathrm{B}}{T}_{\alpha \left(n+1\right)}^i\right]\right)\end{array}$$E83

From Eq. (82) we get

$${\dot{T}}_{\alpha \left(n+1\right)}^i={\mathrm{\gamma }}^{-1}\left[{\dot{T}}_{\mathit{\alpha n}}^i+\frac{\mathrm{\Delta }\tau }{2}\left({\overbrace{\mathrm{X}}}^{-1}\left[\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}^i+\overbrace{\mathbb{R}}-\overbrace{\mathrm{B}}{T}_{\alpha \left(n+1\right)}^i\right]+{\ddot{T}}_{\mathit{\alpha n}}^i\right)\right]$$E84

where $\mathrm{\gamma }=\left(I+\frac{1}{2}\overbrace{\mathrm{A}}\mathrm{\Delta }\tau {\overbrace{\mathrm{X}}}^{-1}\right)$.

Substituting Eq. (84) into Eq. (83), we obtain

$$\begin{array}{ll}{T}_{\alpha \left(n+1\right)}^i=& T_{\mathit{\alpha n}}^i+\mathrm{\Delta }\tau {\dot{T}}_{\mathit{\alpha n}}^i+\frac{\mathrm{\Delta }{\tau }^2}{4}\left({\ddot{T}}_{\mathit{\alpha n}}^i+{\overbrace{\mathrm{X}}}^{-1}\left[\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}^i+\overbrace{\mathbb{R}}-\overbrace{\mathrm{A}}\left({\mathrm{\gamma }}^{-1}\left[{\dot{T}}_{}^{}\right.\right.\right.\right.\\ & \left.\left.\left.\left.+\frac{\mathrm{\Delta }\tau }{2}\left({\overbrace{\mathrm{X}}}^{-1}\left[\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}^i+\overbrace{\mathbb{R}}-\overbrace{\mathrm{B}}{T}_{\alpha \left(n+1\right)}^i\right]+{\ddot{T}}_{\mathit{\alpha n}}^i\right)\right]\right)-\overbrace{\mathrm{B}}{T}_{\alpha \left(n+1\right)}^i\right]\right)\end{array}$$E85

Substituting ${\dot{T}}_{n+1}^i$ from Eq. (84) into Eq. (76), we get

$$\begin{array}{ll}{\ddot{T}}_{\alpha \left(n+1\right)}^i=& {\overbrace{\mathrm{X}}}^{-1}\left[\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}^i+\overbrace{\mathbb{R}}-\overbrace{\mathrm{A}}\left({\mathrm{\gamma }}^{-1}\left[{\dot{T}}_{\mathit{\alpha n}}^i+\frac{\mathrm{\Delta }\tau }{2}\left({\overbrace{\mathrm{X}}}^{-1}\left[\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}^i+\overbrace{\mathbb{R}}\right.\right.\right.\right.\right.\\ & \left.\left.\left.\left.\left.-\overbrace{\mathrm{B}}{T}_{\alpha \left(n+1\right)}^i\right]+{\ddot{T}}_{\alpha \left(n+1\right)}^i\right)\right]\right)-\overbrace{\mathrm{B}}{T}_{\alpha \left(n+1\right)}^i\right]\end{array}$$E86

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: $U_{n+1}^{\mathit{ip}}=U_n^i$. The second step is to substitute for ${\dot{U}}_{n+1}^i$ and ${\ddot{U}}_{n+1}^i$ from Eqs. (77) and (75), respectively, in Eq. (85) and solve the resulted equation for the three-temperature fields. The third step is to correct the displacement using the computed three-temperature fields for the Eq. (80). The fourth step is to compute ${\dot{U}}_{n+1}^i$, ${\ddot{U}}_{n+1}^i$, ${\dot{T}}_{\alpha \left(n+1\right)}^i$, and ${\ddot{T}}_{\alpha \left(n+1\right)}^i$ from Eqs. (79), (81), (82), and (86), respectively.

The continuity conditions for temperature, heat flux, displacement, and traction that have been considered in the current chapter can be expressed as

$${\left.T_{\alpha }^i\left(x,z,\tau \right)\right|}_{x=h^i}={\left.T_{\alpha }^{\left(i+1\right)}\left(x,z,\tau \right)\right|}_{x=h^i}$$E87

$${\left.q^i\left(x,z,\tau \right)\right|}_{x=h^i}={\left.q^{\left(i+1\right)}\left(x,z,\tau \right)\right|}_{x=h^i}$$E88

$${\left.u_f^i\left(x,z,\tau \right)\right|}_{x=h^i}={\left.u_f^{\left(i+1\right)}\left(x,z,\tau \right)\right|}_{x=h^i}$$E89

$${\left.{\overline{t}}_a^i\left(x,z,\tau \right)\right|}_{x=h^i}={\left.{\overline{t}}_a^{\left(i+1\right)}\left(x,z,\tau \right)\right|}_{x=h^i}$$E90

where $n$ is the total number of layers, ${\overline{t}}_a$ are the tractions which is defined by ${\overline{t}}_a={\sigma }_{\mathit{ab}}{n}_b$, and $i=1,2,\dots ,n-1$.

The initial and boundary conditions of the present study are

$$u_f^i\left(x,z,0\right)={\dot{u}}_f^i\left(x,z,0\right)=0\text{for}\left(x,z\right)\in R\cup C$$E91

$$u_f^i\left(x,z,\tau \right)={\mathrm{\Psi }}_f\left(x,z,\tau \right)\text{for}\left(x,z\right)\in C_3$$E92

$${\overline{t}}_a^i\left(x,z,\tau \right)={\mathrm{\Phi }}_f\left(x,z,\tau \right)\text{for}\left(x,z\right)\in C_4,\tau >0$$E93

$$T_{\alpha }^i\left(x,z,0\right)=T_{\alpha }^i\left(x,z,0\right)=0\text{for}\left(x,z\right)\in R\cup C$$E94

$$T_{\alpha }^i\left(x,z,\tau \right)=\overline{f}\left(x,z,\tau \right)\text{for}\left(x,z\right)\in C_1,\tau >0$$E95

$$q^i\left(x,z,\tau \right)=\overline{h}\left(x,z,\tau \right)\text{for}\left(x,z\right)\in C_2,\tau >0$$E96

where ${\mathrm{\Psi }}_f$, ${\mathrm{\Phi }}_f$, f, and $\overline{h}$ are prescribed functions, $C=C_1\cup C_2=C_3\cup C_4$, and $C_1\cap C_2=C_3\cap C_4=\cancel{0}$.

4. 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 ${\alpha }^H$ can be written in terms of the homogenized elasticity matrix D^H and the homogenized vector of stress-temperature coefficients ${\beta }^H$ as follows:

$${\alpha }^H={\left(D^H\right)}^{-1}{\beta }^H$$E97

For the material design, the derivative of the homogenized vector of thermal expansion coefficients can be written as

$$\frac{\partial {\alpha }^H}{\partial X_{\mathit{kl}}^m}={\left(D^H\right)}^{-1}\left(\frac{\partial {\beta }^H}{\partial X_{\mathit{kl}}^m}-\frac{\partial D^H}{\partial X_{\mathit{kl}}^m}{\alpha }^H\right)$$E98

where $\frac{\partial D^H}{\partial X_{\mathit{kl}}^m}$ and $\frac{\partial {\beta }^H}{\partial X_{\mathit{kl}}^m}$ for any lth material phase can be calculated using the adjoint variable method [73] as

$$\frac{\partial D^H}{\partial X_{\mathit{kl}}^m}=\frac{1}{\left|\mathrm{\Omega }\right|}{\int }_Y{\left(I-B^m{U}^m\right)}^T\frac{\partial D^m}{\partial X_{\mathit{kl}}^m}\left(I-B^m{U}^m\right)\mathit{dy}$$E99

and

$$\frac{\partial {\beta }^H}{\partial X_{\mathit{kl}}^m}=\frac{1}{\left|Y\right|}{\int }_Y{\left(I-B^m{U}^m\right)}^T\frac{\partial D^m}{\partial X_{\mathit{kl}}^m}\left({\alpha }^m-B^m{\phi }^m\right)\mathit{dy}+\frac{1}{\left|\mathrm{\Omega }\right|}{\int }_Y{\left(I-B^m{U}^m\right)}^T\frac{\partial {\alpha }^m}{\partial X_{\mathit{kl}}^m}\mathit{dy}$$E100

where $\left|\mathrm{\Omega }\right|$ is the volume of the base cell.

5. 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

$$C_{\mathit{\text{pjkl}}}=\left[\begin{array}{cccccc}430.1& 130.4& 18.2& 0& 0& 201.3\\ 130.4& 116.7& 21.0& 0& 0& 70.1\\ 18.2& 21.0& 73.6& 0& 0& 2.4\\ 0& 0& 0& 19.8& -8.0& 0\\ 0& 0& 0& -8.0& 29.1& 0\\ 201.3& 70.1& 2.4& 0& 0& 147.3\end{array}\right]\mathrm{GPa}$$E101

The mechanical temperature coefficient is

$${\beta }_{\mathit{pj}}=\left[\begin{array}{ccc}1.01& 2.00& 0\\ 2.00& 1.48& 0\\ 0& 0& 7.52\end{array}\right]\cdot {10}^6\frac{\mathrm{N}}{{\mathrm{Km}}^2}$$E102

The thermal conductivity tensor is

$$k_{\mathit{pj}}=\left[\begin{array}{ccc}5.2& 0& 0\\ 0& 7.6& 0\\ 0& 0& 38.3\end{array}\right]\mathrm{W}/\mathrm{Km}$$E103

Mass density $\rho =7820\mathrm{kg}/{\mathrm{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‐\text{axis}$. It was shown from this figure that the fraction order parameter has great effects on the total three-temperature sensitivity.

Figure 2.

Variation of the total 3T sensitivity along x-axis.

Figures 3 and 4 show the variation of the displacement components $u_1$ and $u_2$ along the $x‐\text{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.

Figure 3.

Variation of the displacement u₁ sensitivity along x-axis.

Figure 4.

Variation of the displacement u₂ sensitivity along x-axis.

Figures 5–7 show the variation of the thermal stress components ${\sigma }_{11}$, ${\sigma }_{12}$, and ${\sigma }_{22}$, respectively, along the $x‐\text{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.

Figure 5.

Variation of the thermal stress σ₁₁ sensitivity along x-axis.

Figure 6.

Variation of the thermal stress σ₁₂ sensitivity along x-axis.

Figure 7.

Variation of the thermal stress σ₂₂ sensitivity along x-axis.

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, the displacement $u_1$ and thermal stress ${\sigma }_{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.

Figure 8.

Variation of the displacement u₁ sensitivity along x-axis.

Figure 9.

Variation of the thermal stress σ₁₁ waves along x-axis.

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 $\alpha =0.5$ and shown in Figure 11b for $\alpha =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.

Figure 10.

Design domain of a short cantilever beam.

Figure 11.

The final topology of a short cantilever beam: (a) α = 0.5 and (b) α = 1.0.

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 $\alpha =0.5$ and shown in Figure 13b for $\alpha =1.0$ to show the effect of fractional order parameter on the final topology of the multi-material ISMFGA structure.

Figure 12.

Design domain of a MBB beam.

Figure 13.

The final topology of MBB beam: (a) α = 0.5 and (b) α = 1.0.

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 $\alpha =0.5$ and shown in Figure 15b for $\alpha =1.0$.

Figure 14.

Design domain of a roller-supported beam.

Figure 15.

The final topology of a multi-material roller-supported beam: (a) α = 0.5 and (b) α = 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.

Figure 16.

Design domain of a cantilever beam.

Figure 17.

The final topology of a cantilever beam: (a) MMA and (b) BESO.

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

1. Find X^M

2. That minimize $C^M=\frac{1}{2}{\left(P^M\right)}^T{u}^M=\frac{1}{2}{\left(f^{M,\mathit{ter}}+f^{M,\mathit{mec}}\right)}^T{u}^M$

3. Subject to $V^{M_{,}^{\ast }}-{\sum }_{i=1}^N{V}_i^M{X}_i^M=0$

$$K^M{u}^M=P^M$$

$$X_i^M=x_{\mathit{min}}\mathrm{V}1$$

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_{,}^{\ast }}$ is the volume of the solid material; N is the total number of elements; $K^M$ is the global stiffness matrix; $x_{\mathit{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,\mathit{mec}}$ is the mechanical load vector; and $f^{M,\mathit{ter}}$ is the thermal load vector. Also, the BESO parameters considered in Examples 1–4 can be seen in Tables 1–4, respectively.

Variable name	Variable description	Variable value
$V_f^M$	Final volume fraction	0.5
ER^M	Evolutionary ratio	1%
${\mathit{AR}}_{\mathit{max}}^M$	Volume addition ratio	5%
$r_{\mathit{min}}^M$	Filter ratio	3 mm
τ	Convergence tolerance	0.1%
N	Convergence parameter	5

Table 1.

BESO parameters for minimization of a short cantilever beam.

Variable name	Variable description	Variable value
$V_{f1}^M$	Final volume fraction of the material 1 for both interpolations	0.10
$V_{f2}^M$	Final volume fraction of the material 2 for both interpolations	0.20
ER^M	Evolutionary ratio for interpolation 1	2%
ER^M	Evolutionary ratio for interpolation 2	3%
${\mathit{AR}}_{\mathit{max}}^M$	Volume addition ratio for interpolation 1	3%
${\mathit{AR}}_{\mathit{max}}^M$	Volume addition ratio for interpolation 2	2%
$r_{\mathit{min}}^M$	Filter ratio for interpolation 1	4 mm
$r_{\mathit{min}}^M$	Filter ratio for interpolation 2	3 mm
τ	Convergence tolerance for both interpolations	0.01%
N	Convergence parameter for both interpolations	5

Table 2.

Multi-material BESO parameters for minimization of a MBB beam.

Variable name	Variable description	Variable value
$V_{f1}^M$	Final volume fraction of the material 1 for both interpolations	0.25
$V_{f1}^M$	Final volume fraction of the material 2 for both interpolations	0.25
ER^M	Evolutionary ratio for interpolation 1	3%
ER^M	Evolutionary ratio for interpolation 2	3%
${\mathit{AR}}_{\mathit{max}}^M$	Volume addition ratio for interpolation 1	1%
${\mathit{AR}}_{\mathit{max}}^M$	Volume addition ratio for interpolation 2	1%
$r_{\mathit{min}}^M$	Filter ratio for interpolation 1	4 mm
$r_{\mathit{min}}^M$	Filter ratio for interpolation 2	4 mm
τ	Convergence tolerance for both interpolations	0.5 %
N	Convergence parameter for both interpolations	5

Table 3.

Multi-material BESO parameters for minimization of a roller-supported beam.

Variable name	Variable description	Variable value
$V_f^M$	Final volume fraction	0.4
ER^M	Evolutionary ratio	1.2%
${\mathit{AR}}_{\mathit{max}}^M$	Volume addition ratio	3%
$r_{\mathit{min}}^M$	Filter ratio	0.19 mm
τ	Convergence tolerance	0.1%
N	Convergence parameter	5

Table 4.

BESO parameters for minimization of a cantilever beam.

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.

1. Find X^M

2. That minimize $C^M=\frac{1}{2}{\left(P^M\right)}^T{u}^M=\frac{1}{2}{\left(f^{M,\mathit{ter}}+f^{M,\mathit{mec}}\right)}^T{u}^M$

3. Subject to $V_j^{M_{,}^{\ast }}-{\sum }_{i=1}^N{V}_i^M{X}_{\mathit{ij}}^M-{\sum }_{i=1}^{j-1}{V}_i^{M_{,}^{\ast }}=0;j=1,2$

$$K^M{u}^M=P^M$$

$$X_i^M=x_{\mathit{min}}\mathrm{V}1;j=1,2$$

where $V_j^{M_{,}^{\ast }}$ is the volume of jth material phase and i and j denote the element ith which is made of jth material.

6. 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.