A New Boundary Element Formulation for Modeling and Optimization of Three-Temperature Nonlinear Generalized Magneto-Thermoelastic Problems of FGA Composite Microstructures

Abstract

The main purpose of this chapter is to propose a new boundary element formulation for the modeling and optimization of three-temperature nonlinear generalized magneto-thermoelastic functionally graded anisotropic (FGA) composite microstructures’ problems, which is the gap of this study. Numerical results show that anisotropy and the functionally graded material have great influences on the nonlinear displacement sensitivities and nonlinear thermal stress sensitivities of composite microstructure optimization problem. Since, there are no available data for comparison, except for the problems with one-temperature heat conduction model, we considered the special case of our general study based on replacing three-temperature radiative heat conductions with one-temperature heat conduction. In the considered special case, numerical results demonstrate the validity and accuracy of the proposed technique. In order to solve the optimization problem, the method of moving asymptotes (MMA) based on the bi-evolutionary structural optimization method (BESO) has been implemented. A new class of composite microstructures problems with holes or inclusions was studied. The two-phase magneto-thermoelastic composite microstructure which is studied in this chapter consists of two different FGA materials. Through this chapter, we investigated that the optimal material distribution of the composite microstructures depends strongly on the heat conduction model, functionally graded parameter, and shapes of holes or inclusions.

1. Introduction

In the last few years, there is significant interest in using advanced composite structures, and among the oldest examples of them, reinforced concrete, mixing concrete and steel, and plastics laminated with wood. The main benefit of the composite structures which consist of two or more different materials is that the properties of each material can be combined to form a single unit that performs better than the separate component parts. The most common form of a composite structure in construction is a steel and concrete composite, where concrete works well in pressure but has less resistance to tension. However, steel is extremely strong in tension, and when tied together, it results in a highly efficient and lightweight unit usually used for structures such as buildings and multistory bridges. Although fiberglass and carbon/epoxy composites are not yet as important as the oldest advanced composite structures in terms of tonnage or total revenue, they are very important in engineering, aerospace, transportation, bioengineering, optics, electronics, commodities, chemical plant, and energy industries, especially for the new airplanes that will concentrate on achieving major improvements in the fuel use, emissions, noise, transportation energy consumption, and other important issues to conserve the environment [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21].

Microstructure has been known to play a major role in determining the behavior of material. Therefore, material engineers strive to control the microstructure by improving their properties with the aim of producing a uniform microstructure throughout the material. They also produced FGMs whose microstructures depend on the position by treating the microstructure as a position-dependent variable; the properties of different materials can be combined into one component to achieve an optimum performance in a specific application [22, 23].

In recent years, great attention has been directed toward the study of nonlinear generalized magneto-thermoelastic interactions in functionally graded anisotropic (FGA) structures due to its many applications in physics, geophysics, earthquake engineering, astronautics, aeronautics, mining engineering, military technologies, plasma, robotics, high-energy particle accelerators, nuclear reactors, automobile industries, nuclear plants, soil dynamics, and other engineering and industrial applications. Duhamel [24] and Neuman [25] proposed the classical thermoelasticity (CTE) theory which has the following two paradoxes: first, the infinite propagation speeds of thermal signals are predicted, and second, there is no any elastic term included in heat equation. Biot [26] invented the classical coupled thermoelasticity (CCTE) theory to beat the first paradox in CTE, but CTE and CCTE share the second paradox. Then, numerous generalized thermoelasticity theories have been introduced to overcome the two paradoxes inherent in CTE, such as the extended thermoelasticity (ETE) theory of Lord and Shulman [27]; temperature-rate-dependent thermoelasticity (TRDTE) theory of Green and Lindsay [28]; three linear generalized thermoelasticity theories of Green and Naghdi (GN) [29, 30]; namely I, II, and III, respectively [where, GN theory I is based on Fourier’s law of heat conduction and is identical to CTE theory, GN theory II characterizes the thermoelasticity without energy dissipation (TEWOED), and GN theory III characterizes the thermoelasticity with energy dissipation (TEWED)]; dual phase-lag thermoelasticity (DPLTE) [31, 32]; and three-phase-lag thermoelasticity (TPLTE) [33].

A large amount of research has been done on the generalized problems of thermoelasticity [34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. Our interest in studying the three-temperature thermoelasticity [45, 46, 47, 48, 49] has increased due to its important low-temperature and high-temperature applications. Due to the computational difficulties, inherent in solving three-temperature nonlinear generalized magneto-thermoelastic problems of FGA composite microstructures, the problems become too complicated with no general analytical solution. Therefore, we propose a new boundary element modeling technique which has recently been successfully developed and implemented to obtain the approximate solutions for such problems. Now, the boundary element method (BEM), which is also called boundary integral equation method, has been widely adopted in a large variety of engineering and industrial applications. In the BEM, only the boundary of the solution domain needs to be discretized, so, it has a major advantage over other methods which require the whole domain discretization, such as the finite difference method (FDM) [50, 51, 52], discontinuous Galerkin method (DGM) [53], and finite element method (FEM) [54, 55, 56, 57]. This advantage of BEM over domain methods has significant importance for modeling of nonlinear generalized thermoelastic problems which can be implemented using BEM with little cost and less input data [58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71]. Recently, scientists were convinced that only the FEM method could solve complex engineering problems. But now after the huge achievements of the BEM and its ability to solve complex engineering problems with high efficiency, it gets them to change their conviction. Also, they tried to combine FEM and BEM in the solution of their complex problems.

The main aim of this chapter is to propose a novel boundary element formulation for modeling and optimization of three-temperature nonlinear generalized thermoelastic problems of functionally graded anisotropic (FGA) composite microstructures. The proposed boundary element technique has been implemented successfully for solving several engineering, scientific and industrial applications due to its simplicity, efficiency, ease of use, and applicability [72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85]. The numerical results are presented graphically to show the influence of anisotropy and functionally graded materials on the sensitivities of displacements and thermal stresses. Also, numerical results show the effect of heat conduction model, functionally graded parameter, holes shape, and inclusions shape. Numerical results demonstrate the validity and accuracy of our proposed BEM formulation and technique.

A brief summary of the chapter is as follows: Section 1 introduces an overview of the historical background for a better understanding of the nonlinear generalized magneto-thermoelastic problems and composite materials applications. Section 2 describes the physical modeling of the three-temperature nonlinear generalized thermoelastic problems of FGA composite microstructures. Section 3 outlines the BEM implementation for solving the governing equations of the considered problem to obtain the three temperatures and displacement fields. Section 4 outlines the topology optimization technique used to obtain the optimal composite microstructure with and without holes or inclusions of various shapes. Section 5 presents the new numerical results that describe the effects of anisotropy and functionally graded parameters on the problem’s fields’ sensitivities during the optimization process. Section 6 outlines the significant findings of this chapter.

2. Formulation of the problem

Consider a Cartesian coordinates system ${\mathit{Ox}}_1{x}_2{x}_3$ as shown in Figure 1 . We shall consider a functionally graded anisotropic composite microstructure of a finite thickness $\underset{¯}{\beta }$ placed in a primary magnetic field $H_0$ acting in the direction of the $x_3$ -axis. The considered composite microstructure occupies the region $R=\left\{\left(x_1,x_2,x_3\right):0<x_1<\underset{¯}{\alpha },0<x_2<\underset{¯}{\beta },0<x_3<\underset{¯}{\gamma }\right\}$ with functionally graded material properties in the thickness direction.

Figure 1.

Computational domain of considered structure.

The unified governing equations of three-temperature nonlinear generalized magneto-thermoelasticity for FGA composite microstructures can be expressed as follows [45, 46, 47, 48, 49]:

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

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

$${\tau }_{\mathit{ab}}=\mu {\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)$$ E3

The 2D-3T radiative heat conduction Eqs. (7)–(9) can be expressed as follows:

$$\nabla \left[\left({\delta }_{1j}{\mathbb{K}}_{\alpha }+{\delta }_{2j}{\mathbb{K}}_{\alpha }^{\ast }\right)\nabla T_{\alpha }\left(r,\tau \right)\right]-\overline{\mathbb{W}}\left(r,\tau \right)=C_{\mathit{v\alpha }}\rho {\left(x+1\right)}^m{\delta }_1{\delta }_{1j}\frac{\partial T_{\alpha }\left(r,\tau \right)}{\partial \tau }$$ E4

where

$$\overline{\mathbb{W}}\left(r,\tau \right)=\left\{\begin{array}{lll}\rho {\mathbb{W}}_{\mathit{ei}}\left(T_e-T_i\right)+\rho {\mathbb{W}}_{\mathit{er}}\left(T_e-T_p\right)+\overline{\overline{\mathbb{W}}},\alpha =e,{\delta }_1=1& & \\ -\rho {\mathbb{W}}_{\mathit{ei}}\left(T_e-T_i\right)+\overline{\overline{\mathbb{W}}},& \alpha =i,{\delta }_1=1& ,C_{\mathit{v\alpha }}=\left\{\begin{array}{ll}{c}_e& \alpha =e\\ c_i& \alpha =i\\ c_p{T}_p^3& \alpha =p\end{array}\right.\\ -\rho {\mathbb{W}}_{\mathit{er}}\left(T_e-T_p\right)+\overline{\overline{\mathbb{W}}},& \alpha =p,{\delta }_1=T_p^3& \end{array}\right.$$ E5

in which

$$\begin{array}{l}\overline{\overline{\mathbb{W}}}\left(r,\tau \right)=-{\delta }_{2j}{\mathbb{K}}_{\alpha }{\dot{T}}_{\alpha ,\mathit{ab}}+{\beta }_{\mathit{ab}}{T}_{\alpha 0}{\left(x+1\right)}^m⌈\AA {\delta }_{1i}{\dot{u}}_{a,b}+\left({\tau }_0+{\delta }_{2i}\right){\ddot{u}}_{a,b}⌉\\ +\rho c_{\alpha }{\left(x+1\right)}^m\left[\left({\tau }_0+{\delta }_{1j}{\tau }_2+{\delta }_{2j}\right){\ddot{T}}_{\alpha }\right]\end{array}$$ E6

and

$${\mathbb{W}}_{\mathit{ei}}=\rho {\mathbb{A}}_{\mathit{ei}}{T}_e^{-2/3},{\mathbb{W}}_{\mathit{er}}=\rho {\mathbb{A}}_{\mathit{er}}{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}}$$ E7

where ${\sigma }_{\mathit{ab}}$ , ${\tau }_{\mathit{ab}}$ , $u_k$ , $T_{\alpha }$ , and $T_{\alpha 0}$ are the mechanical stress tensor, Maxwell’s electromagnetic stress tensor, displacement vector, temperature, and reference temperature, respectively; $C_{\mathit{\text{abfg}}}\left(C_{\mathit{\text{abfg}}}=C_{\mathit{\text{fgab}}}=C_{\mathit{\text{bafg}}}\right)$ and ${\beta }_{\mathit{ab}}\left({\beta }_{\mathit{ab}}={\beta }_{\mathit{ba}}\right)$ are, respectively, the constant elastic moduli and stress-temperature coefficients of the anisotropic medium; $\mu$ is the magnetic permeability; $\tilde{h}$ is the perturbed magnetic field; ${\mathbb{K}}_{\alpha }\left(\alpha =e,i,p\right)$ are the thermal conductivity coefficients; $C_{\mathit{v\alpha }}\left(\alpha =e,i,p\right)$ are specific heat coefficients; $e,i$ , and $p$ denote electron, ion, and phonon, respectively; ${\mathbb{K}}_{\alpha }^{\ast }$ is the second order tensor associated with the TEWED and TEWOED theories; ${\mathbb{W}}_{\mathrm{ei}}$ is the electron-ion energy coefficient; ${\mathbb{W}}_{\mathit{ep}}$ is the electron-phonon energy coefficient; $c_{\alpha }\left(\alpha =e,i,p\right)$ are constants; $\rho$ , $\tau$ , and $\AA$ are the density, time, and unified parameter which introduced to consolidate all theories into a unified equations system, respectively; ${\tau }_0$ , ${\tau }_1$ , and ${\tau }_2$ are the relaxation times; and $m$ is a functionally graded parameter. Also, $g_1$ , $g_2$ , ${\mathrm{\Psi }}_f$ , and ${\delta }_f$ are suitably prescribed functions; ${\overline{t}}_a$ are the tractions defined by ${\overline{t}}_a={\sigma }_{\mathit{ab}}{n}_b$ ; and ${\delta }_{1j}$ and ${\delta }_{2j}$ are the Kronecker delta functions.

A superposed dot denotes the differentiation with respect to the time, and a comma followed by a subscript denotes partial differentiation with respect to the corresponding coordinates.

The unit mass total energy can be written as

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

By using the following initial and boundary conditions:

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

$${\left.{\mathbb{K}}_{\alpha }{\left.\frac{\partial T_{\alpha }}{\partial n}\right|}_{C_1}=0,\alpha =e,i,T_r\right|}_{C_1}=g_2\left(x,\tau \right)$$ E10

$${\mathbb{K}}_{\alpha }{\left.\frac{\partial T_{\alpha }}{\partial n}\right|}_{C_2}=0,\alpha =e,i,p$$ E11

$$u_f\left(x,y,0\right)={\dot{u}}_f\left(x,y,0\right)=0\text{for}\left(x,y\right)\in R\cup C$$ E12

$$u_f\left(x,y,\tau \right)={\mathrm{\Psi }}_f\left(x,y,\tau \right)\text{for}\left(x,y\right)\in C_3$$ E13

$${\overline{t}}_a\left(x,y,\tau \right)={\delta }_f\left(x,y,\tau \right)\text{for}\left(x,y\right)\in C_4,\tau >0,C=C_3\cup C_4,C_3\cap C_4=\varnothing$$ E14

By using the fundamental solution that satisfies the following equation:

$$D{\nabla }^2{T}_{\alpha }+\frac{\partial T_{\alpha }^{\ast }}{\partial n}=-\delta \left(r-p_i\right)\delta \left(\tau -r\right),D=\frac{{\mathbb{K}}_{\alpha }}{\mathit{\rho c}}$$ E15

where $p_i$ are singular points.

The above governing Eqs. (1)–(4) can be reduced to the different theories of three-temperature nonlinear generalized magneto-thermoelasticity for FGA composite microstructures as follows [77]:

$$\mathrm{CTE}:j=1,\AA =0\text{and}{\tau }_0={\tau }_1={\tau }_2=0$$ E16

$$\text{CCTE}:j=1,\AA =1\text{and}{\tau }_0={\tau }_1={\tau }_2=0$$ E17

$$\mathrm{ETE}:j=1,\AA =1\text{and}{\tau }_1={\tau }_2=0$$ E18

$$\text{TRDTE}:j=1,\AA =1\text{and}{\tau }_0=0$$ E19

$$\text{TEWED}:j=2,\AA =0\text{and}{\tau }_0=0$$ E20

$$\text{TEWOED}:j=2,\AA =0,{\tau }_0=0\text{and}{\mathbb{K}}_{\alpha }\to 0$$ E21

3. BEM implementation

By using Eqs. (2) and (3), we can write Eq. (1) as follows:

$$L_{\mathit{gb}}{u}_f=\rho {\ddot{u}}_a-D_{a}T=f_{\mathit{gb}}$$ E22

where

$$L_{\mathit{gb}}=D_{\mathit{abf}}\frac{\partial }{\partial x_b}+D_{\mathit{af}}+\mathrm{\Lambda }{D}_{a1f},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}}\left(\frac{\partial }{\partial x_b}+{\delta }_{\mathrm{b}1}\mathrm{\Lambda }+{\tau }_1\left(\frac{\partial }{\partial x_b}+\mathrm{\Lambda }\right)\frac{\partial }{\partial \tau }\right),$$

$$\mathrm{\Lambda }=\frac{m}{x+1},f_{\mathit{gb}}=\rho {\ddot{u}}_a-D_{a}T.$$ E23

The field equations can be written in the following operator form:

$$L_{\mathit{gb}}{u}_f=f_{\mathit{gb}}$$ E24

$$L_{\mathit{ab}}T=f_{\mathit{ab}}$$ E25

where the operators $L_{\mathit{gb}}$ and $f_{\mathit{gb}}$ are defined above in Eq. (23), and the operators $L_{\mathit{ab}}$ and $f_{\mathit{ab}}$ are defined as follows:

$$L_{\mathit{ab}}=\nabla \left({\delta }_{2j}{\mathbb{K}}_{\alpha }^{\ast }\right)\nabla$$ E26

$$f_{\mathit{ab}}=-\nabla \left({\delta }_{1j}{\mathbb{K}}_{\alpha }\right)\nabla +c_{\alpha }{\mathit{\rho \delta }}_1{\delta }_{1j}{\left(x+1\right)}^m\frac{\partial T_{\alpha }\left(r,\tau \right)}{\partial \tau }+\overline{\mathbb{W}}\left(r,\tau \right)$$ E27

By applying the weighted residual method (WRM) to the differential Eq. (24), we obtain

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

Now, we can choose the fundamental solution $u_{\mathit{df}}^{\ast }$ as weighting function as

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

The corresponding traction field can be expressed as

$$t_{\mathit{da}}^{\ast }=C_{\mathit{\text{abfg}}}{\left(x+1\right)}^m{u}_{\mathit{df},g}^{\ast }{n}_b$$ E30

The traction vector can be expressed as

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

By applying integration by parts to Eq. (28) and using the sifting property with Eqs. (29) and (31), we obtain

$$u_d\left(\xi \right)=\underset{C}{\int }\left(u_{\mathit{da}}^{\ast }{t}_a-t_{\mathit{da}}^{\ast }{u}_a+u_{\mathit{da}}^{\ast }{\beta }_{\mathit{ab}}{\left(x+1\right)}^m{\mathit{Tn}}_b\right)\mathit{dC}-\underset{R}{\int }{f}_{\mathit{gb}}{u}_{\mathit{da}}^{\ast }\mathit{dR}$$ E32

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

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

By implementing the WRM and integration by parts, we can write Eq. (25) in the following form:

$$\underset{R}{\int }\left(L_{\mathit{ab}}{\mathit{TT}}^{\ast }-L_{\mathit{ab}}{T}^{\ast }T\right)\mathit{dR}=\underset{C}{\int }\left(q^{\ast }T-{\mathit{qT}}^{\ast }\right)\mathit{dC}$$ E34

where

$$q=-{\mathbb{K}}_{\alpha }{t}_{.b}{n}_a$$ E35

$$q^{\ast }=-{\mathbb{K}}_{\alpha }{T}_{.b}^{\ast }{n}_a$$ E36

Based on the sifting property, we can express Eq. (34) as follows:

$$T\left(\xi \right)=\underset{C}{\int }\left(q^{\ast }T-{\mathit{qT}}^{\ast }\right)\mathit{dC}-\underset{R}{\int }{f}_{\mathit{ab}}{T}^{\ast }\mathit{dR}$$ E37

The field Eqs. (32) and (37) can be written in one equation of the form:

$$\begin{array}{ll}\left[\begin{array}{l}{u}_d\left(\xi \right)\\ T\left(\xi \right)\end{array}\right]= & \underset{C}{\int }\left\{-\left[\begin{array}{cc}{t}_{\mathit{da}}^{\ast }& -u_{\mathit{da}}^{\ast }{\beta }_{\mathit{ab}}{\left(x+1\right)}^m{n}_b\\ 0& -q^{\ast }\end{array}\right]\left[\begin{array}{l}{u}_a\\ T\end{array}\right]+\left[\begin{array}{cc}{u}_{\mathit{da}}^{\ast }& 0\\ 0& -T^{\ast }\end{array}\right]\left[\begin{array}{l}{t}_a\\ q\end{array}\right]\right\}\mathit{dC}\\ & -\underset{R}{\int }\left[\begin{array}{cc}{u}_{\mathit{da}}^{\ast }& 0\\ 0& -T^{\ast }\end{array}\right]\left[\begin{array}{l}{f}_{\mathit{gb}}\\ -f_{\mathit{ab}}\end{array}\right]\mathit{dR}\end{array}$$ E38

The generalized thermoelastic vectors and tensors can be written in contracted notation as follows:

$$U_A=\left\{\begin{array}{ll}{u}_a& a=A=1,2,3\\ T& A=4\end{array}\right.$$ E39

$$T_A=\left\{\begin{array}{ll}{t}_a& a=A=1,2,3\\ q& A=4\end{array}\right.$$ E40

$$U_{\mathit{DA}}^{\ast }=\left\{\begin{array}{ll}{u}_{\mathit{da}}^{\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^{\ast }& D=4;A=4\end{array}\right.$$ E41

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

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

The thermoelastic representation formula (38) can be written in contracted notation as

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

The vector $S_A$ can be splitted as

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

where $S_A^T={\mathrm{\omega }}_{\mathit{AF}}{U}_F$

with

$${\mathrm{\omega }}_{\mathit{AF}}=\left\{\begin{array}{ll}-D_a& A=1,2,3;F=4\\ \nabla \left({\delta }_{2j}{\mathbb{K}}_{\alpha }^{\ast }\right)\nabla T_{\alpha }+\left\{\begin{array}{ll}\rho {\mathbb{W}}_{\mathit{ei}}\left(T_e-T_i\right)+\rho {\mathbb{W}}_{\mathit{er}}\left(T_e-T_p\right),& \alpha =e,{\delta }_1=1\\ -\rho {\mathbb{W}}_{\mathit{ei}}\left(T_e-T_i\right),& \alpha =i,{\delta }_1=1\\ -\rho {\mathbb{W}}_{\mathit{er}}\left(T_e-T_p\right),& \alpha =p,{\delta }_1=\frac{4}{\rho }{T}_p^3\end{array}\right.& \text{otherwise}\end{array}\right.$$ E46

$$S_A^{\dot{T}}=\left(-{\delta }_{2j}{\mathbb{K}}_{\alpha }\frac{\partial }{\partial x_a}\frac{\partial }{\partial x_b}+c_{\alpha }{\mathit{\rho \delta }}_1{\delta }_{1j}{\left(x+1\right)}^m\right){\delta }_{\mathit{AF}}{\dot{U}}_F$$

$$\text{with}{\delta }_{\mathit{AF}}=\left\{\begin{array}{ll}1& A=4;F=4\\ 0& \text{otherwise}\end{array}\right.$$ E47

$$S_A^{\ddot{T}}=-\rho c_{\alpha }{\left(x+1\right)}^m\left({\tau }_0+{\delta }_{1j}{\tau }_2+{\delta }_{2j}\right){\delta }_{\mathit{AF}}{\ddot{U}}_F$$ E48

$$S_A^{\dot{u}}=-{\beta }_{\mathit{ab}}{\left(x+1\right)}^m{T}_{\alpha 0}\AA {\delta }_{1j}{\dot{U}}_F$$ E49

$$S_A^{\ddot{u}}=\mathrm{Ⅎ}{\ddot{U}}_F\text{with}\mathrm{Ⅎ}=\left\{\begin{array}{ll}\rho {\left(x+1\right)}^m& A=1,2,3;\mathrm{F}=1,2,3,\\ -T_{\alpha 0}{\beta }_{\mathit{fg}}{\left(x+1\right)}^m\left({\tau }_0+{\delta }_{2j}\right)A=4;f=F=4& \end{array}\right.$$ E50

The thermoelastic representation formula (38) can also be expressed as follows:

$$\begin{array}{ll}\left[S_A\right]= & \left[\begin{array}{l}-D_a{T}_{\alpha }\\ \nabla \left({\delta }_{2j}{\mathbb{K}}_{\alpha }^{\ast }\right)\nabla T_{\alpha }+\left\{\begin{array}{cc}\rho {\mathbb{W}}_{\mathit{ei}}\left(T_e-T_i\right)+\rho {\mathbb{W}}_{\mathit{er}}\left(T_e-T_p\right),& \alpha =e,{\delta }_1=1\\ -\rho {\mathbb{W}}_{\mathit{ei}}\left(T_e-T_i\right),& \alpha =i,{\delta }_1=1\\ -\rho {\mathbb{W}}_{\mathit{er}}\left(T_e-T_p\right),& \alpha =p,{\delta }_1=\frac{4}{\rho }{T}_p^3\end{array}\right.\end{array}\right]\end{array}$$

$$\begin{array}{l}+\left({\delta }_{2j}{\mathbb{K}}_{\alpha }\frac{\partial }{\partial x_a}\frac{\partial }{\partial x_b}-c_{\alpha }{\mathit{\rho \delta }}_1{\delta }_{1j}{\left(x+1\right)}^m\right)\left[\begin{array}{l}0\\ {\dot{T}}_{\alpha }\end{array}\right]\\ -\rho c_{\alpha }{\left(x+1\right)}^m\left({\tau }_0+{\delta }_{1j}{\tau }_2+{\delta }_{2j}\right)\left[\begin{array}{l}0\\ {\ddot{T}}_{\alpha }\end{array}\right]-{\beta }_{\mathit{ab}}{\left(x+1\right)}^m{T}_{\alpha 0}\AA {\delta }_{1j}\left[\begin{array}{l}0\\ {\dot{u}}_{f,g}\end{array}\right]\\ +\left[\begin{array}{l}\rho {\left(x+1\right)}^m{\ddot{u}}_a\\ -T_{\alpha 0}{\beta }_{\mathit{fg}}{\left(x+1\right)}^m\left({\tau }_0+{\delta }_{2j}\right){\ddot{u}}_{f,g}\end{array}\right]\end{array}$$ E51

In order to transform the domain integral in Eq. (44) to the boundary, we approximate the source vector $S_A$ by a series of given known functions $f_{\mathit{AE}}^q$ and unknown coefficients ${\alpha }_E^q$ :

$$S_A\approx \sum _{q=1}^N{f}_{\mathit{\text{AN}}}^q{\alpha }_N^q$$ E52

Thus, the thermoelastic representation formula (44) can be expressed as

$$U_D\left(\xi \right)=\underset{C}{\int }\left(U_{\mathit{DA}}^{\ast }{T}_A-{\tilde{T}}_{\mathit{DA}}^{\ast }{U}_A\right)\mathit{dC}-\sum _{q=1}^N\underset{R}{\int }{U}_{\mathit{DA}}^{\ast }{f}_{\mathit{\text{AN}}}^q\mathit{dR}{\alpha }_N^q$$ E53

By applying the WRM to the following elastic and thermal equations:

$$L_{\mathit{gb}}{u}_{\mathit{fn}}^q=f_{\mathit{\text{an}}}^q$$ E54

$$L_{\mathit{ab}}{T}^q=f_{\mathit{pj}}^q$$ E55

Now, the weighting functions were chosen as the elastic and thermal fundamental solutions $u_{\mathit{da}}^{\ast }$ and $T^{\ast }$ .

Then, the representation formulae of elastic and thermal fields are given as follows:

$$u_{\mathit{dn}}^q\left(\xi \right)=\underset{C}{\int }\left(u_{\mathit{da}}^{\ast }{t}_{\mathit{\text{an}}}^q-t_{\mathit{da}}^{\ast }{u}_{\mathit{\text{an}}}^q\right)\mathit{dC}-\underset{R}{\int }{u}_{\mathit{da}}^{\ast }{f}_{\mathit{\text{an}}}^q\mathit{dR}$$ E56

$$T^q\left(\xi \right)=\underset{C}{\int }\left(q^{\ast }{T}^q-q^q{T}^{\ast }\right)\mathit{dC}-\underset{R}{\int }{f}^q{T}^{\ast }\mathit{dR}$$ E57

The elastic and thermal representation formulae can be combined in one single equation as

$$U_{\mathit{DN}}^q\left(\xi \right)=\underset{C}{\int }\left(U_{\mathit{DA}}^{\ast }{T}_{\mathit{\text{AN}}}^q-T_{\mathit{DA}}^{\ast }{U}_{\mathit{\text{AN}}}^q\right)\mathit{dC}-\underset{R}{\int }{U}_{\mathit{DA}}^{\ast }{f}_{\mathit{\text{AN}}}^q\mathit{dR}$$ E58

By substituting from Eq. (58) into Eq. (53), the coupled thermoelastic representation formula can be expressed as follows:

$$U_D\left(\xi \right)=\underset{C}{\int }\left(U_{\mathit{DA}}^{\ast }{T}_A-{\check{T}}_{\mathit{DA}}^{\ast }{U}_A\right)\mathit{dC}+\sum _{q=1}^N\left(U_{\mathit{DN}}^q\left(\xi \right)+\underset{C}{\int }\left(T_{\mathit{DA}}^{\ast }{U}_{\mathit{\text{AN}}}^q-U_{\mathit{DA}}^{\ast }{T}_{\mathit{\text{AN}}}^q\right)\mathit{dC}\right){\alpha }_N^q$$ E59

By differentiation of Eq. (59) with respect to ${\xi }_l$ , we obtain

$$\begin{array}{ll}\frac{\partial U_D\left(\xi \right)}{\partial {\xi }_l}=& -\underset{C}{\int }\left(U_{\mathit{DA},l}^{\ast }{T}_A-{\check{T}}_{\mathit{DA},l}^{\ast }{U}_A\right)\mathit{dC}\\ & +\sum _{q=1}^N\left(\frac{\partial U_{\mathit{DN}}^q\left(\xi \right)}{\partial {\xi }_l}-\underset{C}{\int }\left(T_{\mathit{DA},l}^{\ast }{U}_{\mathit{\text{AN}}}^q-U_{\mathit{DA},l}^{\ast }{T}_{\mathit{\text{AN}}}^q\right)\mathit{dC}\right){\alpha }_N^q\end{array}$$ E60

According to the procedure described in Fahmy [78], the boundary integral Eq. (59) can be expressed as

$$\tilde{\zeta }U-\mathit{\eta T}=\left(\zeta \check{U}-\eta \check{\wp }\right)\alpha$$ E61

According to the technique of Partridge et al. [68], the displacements $U_F$ and velocities ${\dot{U}}_F$ can be approximated as

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

$${\dot{U}}_F\approx \sum _{q=1}^N{f}_{\mathit{FD}}^q\left(x\right){\tilde{\gamma }}_D^q$$ E63

where $f_{\mathit{FD}}^q$ are known functions, and ${\gamma }_D^q$ and ${\tilde{\gamma }}_D^q$ are unknown coefficients.

The gradients of the displacement and velocity can be approximated as

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

$${\dot{U}}_{F,g}\approx \sum _{q=1}^N{f}_{\mathit{FD},g}^q\left(x\right){\tilde{\gamma }}_D^q$$ E65

By substituting from Eqs. (62) and (63) into Eqs. (46) and (49), the corresponding source terms can be expressed as

$$S_A^T=\sum _{q=1}^N{S}_{\mathit{AD}}^{T,q}{\gamma }_D^q$$ E66

$$S_A^{\dot{u}}=-{\beta }_{\mathit{ab}}{\left(x+1\right)}^m{T}_{\alpha 0}\AA {\delta }_{1j}\sum _{q=1}^N{S}_{\mathit{AD}}^{\dot{u},q}{\tilde{\gamma }}_D^q$$ E67

where

$$S_{\mathit{AD}}^{T,q}=S_{\mathit{AF}}{f}_{\mathit{FD},q}^q$$ E68

$$S_{\mathit{AD}}^{\dot{u},q}=S_{\mathit{FA}}{f}_{\mathit{FD},g}^q$$ E69

By applying the point collocation procedure of Gaul et al. [10] to Eqs. (52), (62), and (63), we obtain the following equation system:

$$\check{S}=\mathit{J\alpha },U=J^{\prime }\gamma ,\dot{U}=J^{\prime }\tilde{\gamma }$$ E70

Solving the system (70) for $\alpha$ , $\gamma$ , and $\tilde{\gamma }$ yields

$$\alpha =J^{-1}\check{S}\gamma ={J^{\prime }}^{-1}U\tilde{\gamma }={J^{\prime }}^{-1}\dot{U}$$ E71

Now, we can write the coefficients $\alpha$ in terms of nodal values of the displacements, $U$ , velocities, $\dot{U}$ , and accelerations, $\ddot{U}$ as follows:

$$\begin{array}{l}\alpha =J^{-1}\left(\check{S}{}^0+B^T{J}^{\text{'}-1}U\right.\\ +\left[\left({\delta }_{2j}{\mathbb{K}}_{\alpha }\frac{\partial }{\partial x_a}\frac{\partial }{\partial x_b}-c_{\alpha }{\mathit{\rho \delta }}_1{\delta }_{1j}{\left(x+1\right)}^m\right){\delta }_{\mathit{AF}}-{\beta }_{\mathit{ab}}(x\right.\\ \left.+1){}^m{T}_{\alpha 0}\AA {\delta }_{1j}{J}^{\text{'}-1}\right]\dot{U}\\ \left.+\left[-c_{\alpha }\rho {\left(x+1\right)}^m\left({\tau }_0+{\delta }_{1j}{\tau }_2+{\delta }_{2j}\right){\delta }_{\mathit{AF}}\right]\ddot{U}\right)\end{array}$$ E72

By substituting from Eq. (72) into Eq. (61) and implementing implicit-implicit staggered algorithm of Farhat et al. [86], the governing equations can be rewritten as

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

$$\overbrace{\mathrm{X}}\ddot{T}+\overbrace{\mathrm{A}}\dot{T}+\overbrace{\mathrm{B}}T=\overbrace{\mathbb{Z}}\ddot{U}+\overbrace{\mathbb{R}}\dot{U}$$ E74

where

$$V=\left(\eta \check{\wp }-\zeta \check{U}\right)J^{-1},\overbrace{M}=V\tilde{A},\overbrace{\mathrm{X}}=-\rho c_{\alpha }{\left(x+1\right)}^m\left({\tau }_0+{\delta }_{1j}{\tau }_2+{\delta }_{2j}\right),$$

$$\overbrace{K}=\tilde{\zeta }+{\mathit{VB}}^T{J^{\prime }}^{-1},\mathbb{Q}=\mathit{\eta T}+V{\stackrel{⌣}{S}}^0,\overbrace{\mathrm{B}}={\delta }_{1j}{\mathbb{K}}_{\alpha }+{\delta }_{2j}{\mathbb{K}}_{\alpha }^{\ast },$$

$$\overbrace{\mathrm{\Gamma }}=V\left[\left({\mathbb{K}}_{\alpha }\frac{\partial }{\partial x_a}\frac{\partial }{\partial x_b}-c_{\alpha }\rho {\left(x+1\right)}^m{\delta }_{1j}\right){\delta }_{\mathit{AF}}-T_0\AA {\delta }_{1j}{\beta }_{\mathit{fg}}{\left(x+1\right)}^m{J^{\prime }}^{-1}\right],$$

$$\overbrace{\mathbb{R}}=T_0{\beta }_{\mathit{ab}}{\left(x+1\right)}^m\AA {\delta }_{1j},\overbrace{\mathbb{Z}}=T_{\alpha 0}{\beta }_{\mathit{ab}}{\left(x+1\right)}^m\left({\tau }_0+{\delta }_{2j}\right),$$

$$\overbrace{\mathrm{A}}=\left({\delta }_{2j}{\mathbb{K}}_{\alpha }^{\ast }\frac{\partial }{\partial x_a}\frac{\partial }{\partial x_b}-\rho c_{\alpha }{\left(x+1\right)}^m{\delta }_{1j}\right){\delta }_{\mathit{AF}}.$$ E75

where $V$ , $\overbrace{M}$ , $\overbrace{\mathrm{\Gamma }}$ , $\overbrace{K}$ , $\overbrace{\mathrm{A}}$ , and $\overbrace{\mathrm{B}}$ are represent the volume, mass, damping, stiffness, capacity, and conductivity matrices, respectively; $\ddot{U}$ , $\dot{U}$ , $U$ , $T$ , and $\overbrace{\mathbb{Q}}$ represent the acceleration, velocity, displacement, temperature, and external force vectors, respectively, $\overbrace{\mathrm{X}}$ is a Green and Lindsay material constants vector, and $\overbrace{\mathbb{Z}}$ and $\overbrace{\mathbb{R}}$ are coupling matrices.

Hence, the governing equations lead to the following coupled system of differential-algebraic equations (DAEs) as in Farhat et al. [86]:

$$\overbrace{M}{\ddot{U}}_{n+1}+\overbrace{\mathrm{\Gamma }}{\dot{U}}_{n+1}+\overbrace{K}{U}_{n+1}={\overbrace{\mathbb{Q}}}_{n+1}^p$$ E76

$$\overbrace{X}{\ddot{T}}_{n+1}+\overbrace{\mathrm{A}}{\dot{T}}_{n+1}+\overbrace{\mathrm{B}}{T}_{n+1}=\overbrace{\mathbb{Z}}{\ddot{U}}_{n+1}\overbrace{\mathbb{R}}{\dot{U}}_{n+1}$$ E77

where ${\overbrace{\mathbb{Q}}}_{n+1}^p=\eta T_{n+1}^p+V{\check{S}}^0$ and $T_{n+1}^p$ .

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

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

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

From Eq. (78) we obtain

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

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

Substitution of Eq. (80) in Eq. (79), we obtain

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

Substituting ${\dot{U}}_{n+1}$ from Eq. (80) into Eq. (76), we obtain

$${\ddot{U}}_{n+1}={\overbrace{M}}^{-1}⌈{\overbrace{\mathbb{Q}}}_{n+1}^p-\overbrace{\mathrm{\Gamma }}⌈{\overline{\mathrm{\gamma }}}^{-1}\left[{\dot{U}}_n+\frac{\mathrm{\Delta }\tau }{2}\left[{\ddot{U}}_n+{\overbrace{M}}^{-1}\left({\overbrace{\mathbb{Q}}}_{n+1}^p-\overbrace{K}{U}_{n+1}\right)\right]\right]⌉-\overbrace{K}{U}_{n+1}⌉$$ E82

Integrating the heat Eq. (74) using the trapezoidal rule and Eq. (77), we get

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

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

From Eq. (83), we have

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

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

On substitution of Eq. (85) in Eq. (84), we obtain

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

On substitution of ${\dot{T}}_{n+1}$ from Eq. (85) in Eq. (77), we get

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

Now, our algorithm for the solution of Eqs. (81) and (86) is obtained as follows:

First step. Predict the displacement field: $U_{n+1}^p=U_n$ .

Second step. Substituting for ${\dot{U}}_{n+1}$ from Eq. (78) and substituting for ${\ddot{U}}_{n+1}$ from Eq. (76). Then, by using the resulted equations in Eq. (86) to obtain the temperature field.

Third step. Correct the displacement field (81) by using the computed temperature.

Fourth step. Compute ${\dot{U}}_{n+1}$ , ${\ddot{U}}_{n+1}$ , ${\dot{T}}_{n+1}$ , and ${\ddot{T}}_{n+1}$ from Eqs. (80), (82), (85), and (87), respectively.

4. Design sensitivity and optimization

According to Fahmy [77, 78], the design sensitivities of the nonlinear temperature field and nonlinear displacement field can be performed by the implicit differentiation of Eqs. (76) and (77), respectively, which describe the structural response with respect to the design variables, then we can compute the nonlinear thermal stresses sensitivities.

In order to solve our topology optimization problem, the method of moving asymptotes (MMA) [87] has been implemented as an optimizer in our topology optimization program. The benefit of MMA algorithm is that it replaces the original nonlinear, non-convex optimization problem by a sequence of approximating convex subproblems which are much easier to solve. The implemented MMA is based on the bi-directional evolutionary structural optimization (BESO), which is the evolutionary topology optimization approach that allows modification of the structure by either adding efficient material or removing inefficient material to or from the structure design [88, 89, 90, 91, 92, 93, 94, 95, 96]. This addition or removal depends upon the sensitivity analysis. Sensitivity analysis is the estimation of the response of the structure to the modification of the input design variables and is dependent upon the calculation of derivatives.

The homogenized vector of thermal expansion coefficients ${\alpha }^H$ can be written in terms of the homogenized elastic matrix $D^H$ and homogenized stress-temperature coefficients vector ${\beta }^H$ as follows:

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

For the material design, the derivative of the homogenized thermal expansion coefficients vector can be expressed 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)$$ E89

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

$$\frac{\partial D^H}{\partial X_{\mathit{kl}}^m}=\frac{1}{\mid \mathrm{\Omega }\mid }{\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}$$ E90

and

$$\begin{array}{ll}\frac{\partial {\beta }^H}{\partial X_{\mathit{kl}}^m}= & \frac{1}{\mid \mathrm{\Omega }\mid }{\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}{\mid \mathrm{\Omega }\mid }{\int }_Y{\left(I-B^m{U}^m\right)}^T{D}^m\frac{\partial {\alpha }^m}{\partial X_{\mathit{kl}}^m}\mathit{dy}\end{array}$$ E91

where, $\mid \mathrm{\Omega }\mid$ is the volume of the base cell.

5. Numerical examples, results, and discussion

The proposed technique used in the current chapter should be applicable to any three-temperature nonlinear generalized magneto-thermoelastic problem. The application is for the purpose of illustration.

The two anisotropic materials considered in the calculation are monoclinic graphite-epoxy and North Sea sandstone reservoir rock, where the physical data of monoclinic graphite-epoxy material is given as follows:

Elasticity tensor:

$$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}$$ E92

Mechanical temperature coefficient:

$${\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\mathrm{N}/{\mathrm{Km}}^2$$ E93

Tensor of thermal conductivity:

$$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}$$ E94

Mass density $\rho =7820\mathrm{kg}/{\mathrm{m}}^3$ and heat capacity $c=461$ J/(kg⋅K), $H_0=1000000$ Oersted, $\mu =0.5$ Gauss/Oersted, $h=2$ , and $\mathrm{\Delta }\tau =0.0001$ .

The physical data of the North Sea sandstone reservoir rock is given as follows:

Elasticity tensor:

$$C_{\mathit{\text{pjkl}}}=\left[\begin{array}{cc}\begin{array}{ccc}17.77& 3.78& 3.76\\ 3.78& 19.45& 4.13\\ 3.76& 4.13& 21.79\end{array}& \begin{array}{ccc}0.24& -0.28& 0.03\\ 0& 0& 1.13\\ 0& 0& 0.38\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0.03& 1.13& 0.38\end{array}& \begin{array}{ccc}8.30& 0.66& 0\\ 0.66& 7.62& 0\\ 0& 0& 7.77\end{array}\end{array}\right]\mathrm{GPa}$$ E95

Mechanical temperature coefficient:

$${\beta }_{\mathit{pj}}=\left[\begin{array}{ccc}0.001& 0.02& 0\\ 0.02& 0.006& 0\\ 0& 0& 0.05\end{array}\right]\cdot {10}^6\mathrm{N}/{\mathrm{Km}}^2$$ E96

Tensor of thermal conductivity:

$$k_{\mathit{pj}}=\left[\begin{array}{ccc}1& 0.1& 0.2\\ 0.1& 1.1& 0.15\\ 0.2& 0.15& 0.9\end{array}\right]\mathrm{W}/\mathrm{km}$$ E97

Mass density $\rho =2216\mathrm{kg}/{\mathrm{m}}^3$ and heat capacity $c=0.1$ J/(kg⋅K), $H_0=1000000$ Oersted, $\mu =0.5$ Gauss/Oersted, $h=2$ , and $\mathrm{\Delta }\tau =0.0001$ .

The initial and boundary conditions considered in the calculations are

$$\mathrm{at}\tau =0u_1=u_2={\dot{u}}_1={\dot{u}}_2=0,T=0$$ E98

$$\mathrm{at}x=0\frac{\partial u_1}{\partial x}=\frac{\partial u_2}{\partial x}=0,\frac{\partial T}{\partial x}=0$$ E99

$$\mathrm{at}x=h\frac{\partial u_1}{\partial x}=\frac{\partial u_2}{\partial x}=0,\frac{\partial T}{\partial x}=0$$ E100

$$\mathrm{at}y=0\frac{\partial u_1}{\partial y}=\frac{\partial u_2}{\partial y}=0,\frac{\partial T}{\partial y}=0$$ E101

$$\mathrm{at}y=b\frac{\partial u_1}{\partial y}=\frac{\partial u_2}{\partial y}=0,\frac{\partial T}{\partial y}=0$$ E102

In order to study the effects of anisotropy and functionally graded materials on composite microstructure, we consider the following four cases, namely, isotropic homogeneous (IH), isotropic functionally graded (IF), anisotropic homogeneous (AH), and anisotropic functionally graded (AF). Also, we considered total temperature $T\left(T=T_e+T_i+T_p\right)$ as the considered temperature field in all calculations of this study.

Figure 2 shows the variations of the nonlinear three-temperature ${\mathrm{T}}_{\mathrm{e}}$ , ${\mathrm{T}}_{\mathrm{i}}$ , and ${\mathrm{T}}_{\mathrm{p}}$ and total temperature $\mathrm{T}\left(\mathrm{T}={\mathrm{T}}_{\mathrm{e}}+{\mathrm{T}}_{\mathrm{i}}+{\mathrm{T}}_{\mathrm{p}}\right)$ , with the time $\mathrm{\tau }$ through composite microstructure.

Figure 2.

Variation of the temperature sensitivity with time τ.

Figures 3 and 4 show the variation of the nonlinear displacement sensitivities ${\mathrm{u}}_1$ and ${\mathrm{u}}_2$ , with time $\tau$ for different cases IH, IF, AH, and AF. It was shown from these figures that the anisotropy and functionally graded material have great effects on the nonlinear displacement sensitivities through the FGA composite microstructure.

Figure 3.

Variation of the displacement u₁ sensitivity with time τ.

Figure 4.

Variation of the displacement u₂ sensitivity with time τ.

Figures 5 – 7 show the variation of the nonlinear thermal stress sensitivities ${\mathrm{\sigma }}_{11}$ , ${\mathrm{\sigma }}_{12}$ , and ${\mathrm{\sigma }}_{22}$ , respectively, with time $\mathrm{\tau }$ for different cases IH, IF, AH, and AF. It was noted from these figures that the anisotropy and functionally graded material have great influences on the nonlinear thermal stress sensitivities through the FGA composite microstructure.

Figure 5.

Variation of the thermal stress σ₁₁ sensitivity with time τ.

Figure 6.

Variation of the thermal stress σ₁₂ sensitivity with time τ.

Figure 7.

Variation of the thermal stress σ₂₂ sensitivity with time τ.

For comparison purposes with those of other studies, we only considered one-dimensional numerical results of the considered three-temperature problem. In the considered special case, the nonlinear displacement $u_1$ and nonlinear thermal stress ${\sigma }_{11}$ results are plotted in Figures 8 and 9 , respectively. It can be noticed from these that the BEM results, which are based on replacing one-temperature heat conduction with three-temperature heat conduction, are in excellent agreement when compared to results obtained from the finite difference method of Pazera and Jędrysiak [97] and the finite element method (FEM) of Xiong and Tian [98]. We thus demonstrate the validity and accuracy of our proposed BEM technique.

Figure 8.

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

Figure 9.

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

Three numerical examples of BESO topological optimization of composite microstructures are performed to illustrate the optimization results of this study [99]. In order to obtain the functionally graded parameter effects during the optimization process of the considered composite microstructure, we consider the following values $m=0,0.5,0.75$ , and $1$ in the one-temperature heat conduction model and the three-temperature radiative heat conduction model.

Example 1. Composite microstructures without holes or inclusions.

The mean compliance has been minimized, to obtain the maximum stiffness for the composite microstructures made from two competitive materials and without holes or inclusions. Investigation of the effect of the functionally graded parameter on the optimal composite microstructure has been shown in Table 1 for the $1T$ model and in Table 2 for the $3T$ model. It is noticed from these tables that the heat conduction model and functionally graded parameter have a significant effect on the topology optimization process of the multi-material FGA composite microstructures.

Table 1.

Investigation of the influence of functionally graded parameter, m, on the optimal composite microstructure for the 1T model.

Table 2.

Investigation of the influence of functionally graded parameter, m, on the optimal composite microstructure for the 3T model.

Example 2. Composite microstructures with circular or square holes.

The mean compliance has been minimized to obtain the maximum stiffness for the composite microstructures made from two competitive materials and with circular or square holes. Investigation of the effect of the functionally graded parameter on the optimal composite microstructure with circular holes has been shown in Table 3 for the $1T$ model and in Table 4 for the $3T$ model. Also, the investigation of the effect of the functionally graded parameter on the optimal composite microstructure with square holes has been shown in Table 5 for the $1T$ model and in Table 6 for the $3T$ model. It is noticed from these tables that the heat conduction model, functionally graded parameter, and holes shape have a significant effect on the topology optimization process of the multi-material FGA composite microstructures.

Table 3.

Investigation of the influence of functionally graded parameter, m, on the optimal composite microstructure with circular shape holes for the 1T model.

Table 4.

Investigation of the influence of functionally graded parameter m on the optimal composite microstructure with circular shape holes for the 3T model.

Table 5.

Investigation of the influence of functionally graded parameter, m, on the optimal composite microstructure with square shape holes for the 1T model.

Table 6.

Investigation of the influence of functionally graded parameter, m, on the optimal composite microstructure with square shape holes for the 3T model.

Example 3. Composite microstructures with circular or square inclusions.

The mean compliance has been minimized to obtain the maximum stiffness for the composite microstructures made from two competitive materials and with circular or square inclusions. Investigation of the effect of the functionally graded parameter on optimal composite microstructure with circular inclusions has been shown in Table 7 for the $1T$ model and in Table 8 for the $3T$ model. Also, the investigation of the effect of the functionally graded parameter on the optimal composite microstructure with square inclusions has been shown in Table 9 for the $1T$ model and in Table 10 for the $3T$ model. It is noticed from these tables that the heat conduction model, functionally graded parameter, and inclusions shape have a significant effect on the topology optimization process of the multi-material FGA composite microstructures.

Table 7.

Investigation of the influence of functionally graded parameter, m, on the optimal composite microstructure with circular shape inclusions for the 1T model.

Table 8.

Investigation of the influence of functionally graded parameter, m, on the optimal composite microstructure with circular shape inclusions for the 3T model.

Table 9.

Investigation of the influence of the functionally graded parameter, m, on the optimal composite microstructure with square shape inclusions for the 1T model.

Table 10.

Investigation of the influence of functionally graded parameter, m, on the optimal composite microstructure with square shape inclusions for the 3T model.

The BESO topology optimization problem implemented in the numerical examples to find the distribution of the two materials in the design domain that minimize the compliance of the structure subject to a volume constraint in both phases can be stated as

Find $X^M$

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$

Subject to $V_j^{M,\ast }-{\mathrm{\Sigma }}_{i=1}^N{V}_i^M{X}_{\mathit{ij}}^M-{\mathrm{\Sigma }}_{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 $X^M$ is the design variable; $V_j^{M,\ast }$ is the volume of the $j\mathrm{th}$ material phase, where $i$ and $j$ denote the element $i\mathrm{th}$ which is made of $j\mathrm{th}$ material; $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 guarantees 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 these examples can be seen in Table 11 . The validity of our implemented BESO topology optimization technique has been demonstrated in our recent reference [100].

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

Multi-material BESO parameters for minimization of a composite microstructure.

Example 4. Laminated composite microstructure with three different sets of boundary conditions are considered in this example to validate the BEM formulation of the current study. These boundary conditions are called: simply—simply supported (SS), clamped—clamped (CC), and clamped—simply supported (CS). One-temperature (1T) and three-temperature (3T) models of nonlinear thermal stresses sensitivities results have been compared with the finite element method (FEM) results of Rajanna et al. [101] as well as with the finite volume method (FVM) results of Fallah and Delzendeh [102], which are tabulated in Table 12 for different types of boundary conditions and different methods. It can be observed that the BEM results for all the three types of boundary conditions are in excellent agreement with FEM results of [101] and the FVM results of [102].

Model	Type	Method	σ 11 sensitivity	σ 12 sensitivity	σ 22 sensitivity
		BEM (present)	0.4084297	0.0509346	0.5332620
IT	SS	FEM [101]	0.4084297	0.0509346	0.5332620
		FVM [102]	0.4084297	0.0509346	0.5332620
		BEM (present)	0.3591487	0.0408259	0.3758618
IT	CC	FEM [101]	0.3591487	0.0408259	0.3758618
		FVM [102]	0.3591487	0.0408259	0.3758618
		BEM (present)	0.2518379	0.0307736	0.2613532
IT	CS	FEM [101]	0.2518378	0.0307735	0.2613531
		FVM [102]	0.2518379	0.0307736	0.2613532
		BEM (present)	0.3147697	0.0304365	0.4767924
3T	SS	FEM [101]	0.3147696	0.0304364	0.4767923
		FVM [102]	0.3147697	0.0304365	0.4767924
		BEM (present)	0.2432756	0.0204748	0.3052857
3T	CC	FEM [101]	0.2432755	0.0204747	0.3052856
		FVM [102]	0.2432756	0.0204748	0.3052857
		BEM (present)	0.1258948	0.0107825	0.2079735
3T	CS	FEM [101]	0.1258947	0.0107824	0.2079734
		FVM [102]	0.1258948	0.0107825	0.2079737

Table 12.

Models of 1T and 3T nonlinear thermal stresses’ sensitivities for different types of boundary conditions and different methods.

6. Conclusion

The main aim of this chapter is to describe a new boundary element formulation for the modeling and optimization of the three-temperature nonlinear generalized magneto-thermoelastic functionally graded anisotropic (FGA) composite microstructures. The governing equations of the considered model are very difficult to solve analytically because of the nonlinearity and anisotropy. To overcome this, we propose a new boundary element formulation for solving such equations, where we used the three-temperature nonlinear radiative heat conduction equations combined with electron, ion, and phonon temperatures. Numerical results show the three-temperature distributions through composite microstructure. The effects of anisotropy and functionally graded material on the three-temperature nonlinear displacement sensitivities and nonlinear thermal stress sensitivities through the composite microstructure are very significant and pronounced. Because there are no available results in the literature to confirm the validity and accuracy of our proposed technique except for one-temperature heat conduction, we replace the three-temperature radiative heat conduction with one-temperature heat conduction as a special case from our current general study. In the considered special case, the BEM results have been compared graphically with the FDM results and FEM results, and it can be noticed that the BEM results are in excellent agreement with the FDM and FEM results. These results thus demonstrate the validity and accuracy of our proposed technique.

Numerical examples are solved using the method of moving asymptotes (MMA) algorithm based on the bi-evolutionary structural optimization method (BESO), where we used the topological optimization to manufacture three-temperature magneto-thermoelastic composite microstructures to obtain the required specific engineering properties. A new class of FGA composite microstructures consisting of two competitive materials has been studied, taking into account existence of holes or inclusions. The effects of the heat conduction model, functionally graded parameter, and holes shape and inclusions shape on the optimal composite microstructure are investigated through the considered examples with great practical interest.

The ability to understand and manipulate composite microstructures has been fundamental to our technical development over time. Today, scientists and engineers recognize the importance of composite microstructures use for economic and environmental reasons. Based on the BEM implementation and its results, this study concluded that the boundary element technique is the most suitable technique for the manufacturing of FGA composite microstructures in the future works. This technique aimed to describe the behavior of FGA composite microstructures and achieves improvement in the composition optimization and mechanical properties of the resulting FGA composite microstructures.

Due to three-temperature and numerous low-temperature and high-temperature applications in laminated composites microstructures, as a future work and based on the findings obtained in the present study, we would suggest further research to develop numerical techniques for solving the three-temperature nonlinear thermoelastic wave propagation problems and for manufacturing of advanced laminated composites. The numerical results of our considered study can provide data references for mechanical engineers, computer engineers, geotechnical engineers, geothermal engineers, technologists, new materials designers, physicists, material science researchers, and those who are interested in novel technologies in the area of three-temperature magneto-thermoelastic FGA composite microstructures. Application of three-temperature theories in advanced manufacturing technologies, with the development of soft machines and robotics in biomedical engineering and advanced manufacturing, and nonlinear generalized magneto-thermoelastic problems will be encountered more often where three-temperature radiative heat conduction will turn out to be the best choice for thermomechanical analysis in the design and analysis of advanced composite microstructures.