A Novel MDD-Based BEM Model for Transient 3T Nonlinear Thermal Stresses in FGA Smart Structures

Abstract

The main objective of this chapter is to introduce a novel memory-dependent derivative (MDD) model based on the boundary element method (BEM) for solving transient three-temperature (3T) nonlinear thermal stress problems in functionally graded anisotropic (FGA) smart structures. The governing equations of the considered study are nonlinear and very difficult if not impossible to solve analytically. Therefore, we develop a new boundary element scheme for solving such equations. The numerical results are presented highlighting the effects of the MDD on the temperatures and nonlinear thermal stress distributions and also the effect of anisotropy on the nonlinear thermal stress distributions in FGA smart structures. The numerical results also verify the validity and accuracy of the proposed methodology. The computing performance of the proposed model has been performed using communication-avoiding Arnoldi procedure. We can conclude that the results of this chapter contribute to increase our understanding on the FGA smart structures. Consequently, the results also contribute to the further development of technological and industrial applications of FGA smart structures of various characteristics.

1. Introduction

Smart materials, which are also called intelligent materials, are engineered materials that have the ability to respond to the changes that occur around them in a controlled fashion by external stimuli, such as stress, heat, light, ultraviolet, moisture, chemical compounds, mechanical strength, and electric and magnetic fields. We can simply define smart materials as materials which adapt themselves as per required condition. The history of the discovery of these materials dates back to the 1880s when Jacques and Pierre Curie noticed a phenomenon that pressure generates electrification around a number of minerals such as quartz and tourmaline, and this phenomenon is called piezoelectric effect, so the piezoelectric materials are the oldest type of smart materials, which are utilized extensively in the fabrication of various devices such as transducers, sensors, actuators, surface acoustic wave devices, frequency control, etc. There are a lot of smart material types like piezoelectric materials, thermochromic pigments, shape memory alloys, magnetostrictive, shape memory polymers, hydrogels, electroactive polymers and bi-component fibers, etc.

Anisotropic smart structures (ASSs) are getting great attention of researchers due to their applications in textile, aerospace, mass transit, marine, automotive, computers and other electronic industries, consumer goods applications, mechanical and civil engineering, infertility treatment, micropumps, medical equipment applications, ultrasonic micromotors, microvalves and photovoltaics, rotating machinery applications, and much more [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12].

The classical thermoelasticity (CTE) theory of Duhamel [13] and Neumann [14] has two shortcomings based on parabolic heat conduction equation of this theory: the first does not involve any elastic terms, while the second has infinite propagation speeds of thermoelastic waves. In order to overcome the first shortcoming, Biot [15] proposed the classical coupled thermoelasticity (CCTE). But CTE and CCTE have the second shortcoming. So, several generalized thermoelasticity theories have been developed to overcome the second shortcoming of CTE. Among of these theories are Lord and Shulman (LS) [16], Green and Lindsay (GL) [17], and Green and Naghdi [18, 19] theories of thermoelasticity with and without energy dissipation, dual-phase-lag thermoelasticity (DPLTE) [20, 21] and three-phase-lag thermoelasticity (TPLTE) [22]. Although thermoelastic phenomena in the majority of practical applications are adequately modeled with the classical Fourier heat conduction equation, there are an important number of problems that require consideration of nonlinear heat conduction equation. It is appropriate in these cases to apply the nonlinear generalized theory of thermoelasticity; great attention has been paid to investigate the nonlinear generalized thermoelastic problems by using numerical methods [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. Fahmy [35, 36, 37, 38, 39] introduced the mathematical foundations of three-temperature (3T) field to thermoelasticity.

The fractional calculus is the mathematical branch that used to study the theory and applications of derivatives and integrals of arbitrary non-integer order. This branch has emerged in recent years as an effective tool for modeling and simulation of various engineering and industrial applications [40, 41]. Due to the nonlocal nature of fractional order operators, they are useful for describing the memory and hereditary properties of various materials and processes. Also, the fractional calculus has drawn wide attention from the researchers of various countries in recent years due to its applications in solid mechanics, fluid dynamics, viscoelasticity, heat conduction modeling and identification, biology, food engineering, econophysics, biophysics, biochemistry, electrochemistry, electrical engineering, finance and control theory, robotics and control theory, signal and image processing, electronics, electric circuits, wave propagation, nanotechnology, etc. [42, 43, 44].

Several mathematics researchers 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 presented the first formula for the fractional order derivative appeared in 1819, where he introduced the $\mathrm{n}$th derivative of the function $y=x^m$ as follows:

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

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

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

By using Cauchy’s integral formula for complex valued analytical functions, Laurent defined the integration of arbitrary order $v>0$ as

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

where ${}_c{D}_x^v$ denotes differentiation of order $v$ of the function f along the x-axis.

Cauchy presented the following fractional order derivative:

$$f_{+}^{\left(\alpha \right)}=\int f\left(\tau \right)\frac{{\left(t-\tau \right)}^{-\alpha -1}}{\mathrm{\Gamma }\left(-a\right)}\mathit{d\tau }$$E4

In 1967, the Italian mathematician Caputo presented his fractional derivative of order $\alpha >0$ as

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

Diethelm [45] has suggested the Caputo derivative to be in the following form:

$$D_a^{\zeta }f\left(\tau \right)={\int }_a^{\tau }{K}_{\zeta }\left(\tau ‐\xi \right)f^{\left(m\right)}\left(\xi \right)\mathit{d\xi }$$E6

where f^(m) is the mth order derivative and m is an integer such that $m‐1<\zeta \le m$

$$K_{\zeta }\left(\tau ‐\xi \right)=\frac{{\left(\tau ‐\xi \right)}^{m‐\zeta ‐1}}{\mathrm{\Gamma }\left(m‐\zeta \right)}$$E7

Wang and Li [46] have introduced a memory-dependent derivative (MDD)

$$D_{\omega }^{\zeta }f\left(\tau \right)=\frac{1}{\omega }{\int }_{\tau ‐\omega }^{\tau }{K}_{\zeta }\left(\tau ‐\xi \right)f^{\left(m\right)}\left(\xi \right)\mathit{d\xi }$$E8

where the first-order $\left(\zeta =1\right)$ of MDD for a differentiable function $f\left(\tau \right)$ can be expressed as

$$D_{\omega }f\left(\tau \right)=\frac{1}{\omega }{\int }_{\tau ‐\omega }^{\tau }K\left(\tau ‐\xi \right)f^{\prime }\left(\xi \right)\mathit{d\xi }$$E9

Based on several practical applications, the memory effect needs weight $0\le K\left(\tau ‐\xi \right)<1$ for $\xi \in \left[\tau ‐\omega ,\tau \right]$, so the MDD magnitude $D_{\omega }f\left(\tau \right)$ is usually smaller than $f^{\prime }\left(\tau \right)$, where the time delay $\left(\omega >0\right)$ and the kernel function ($0\le K\left(\tau ‐\xi \right)\le 1$ for $\xi \in \left[\tau ‐\xi ,\tau \right]$) can be chosen arbitrarily on the delayed interval $\left[\tau ‐\omega ,\tau \right]$, the practical kernel functions are $1,\left[1-\left(\tau ‐\xi \right)\right]$ and ${\left[1-\frac{\tau ‐\xi }{\omega }\right]}^p$, $p=\frac{1}{4},1,2$, etc. These functions are monotonically increasing with $K=0$ for the past time $\tau ‐\xi$ and $K=1$ for the present time $\mathrm{\tau }$. The main feature of MDD is that the real-time functional value depends also on the past time $[\tau ‐\xi ‐\tau [$. So, $D_{\omega }$ depends on the past time (nonlocal operator), while the integration does not depend on the past time (local operator).

As a special case $K\left(\tau ‐\xi \right)\equiv 1$, we have

$${\mathrm{D}}_{\mathrm{\omega }}\mathrm{f}\left(\mathrm{\tau }\right)=\frac{1}{\mathrm{\omega }}{\int }_{\mathrm{\tau }-\mathrm{\omega }}^{\mathrm{\tau }}{\mathrm{f}}^{\prime }\left(\mathrm{\xi }\right)\mathrm{d\xi }=\frac{\mathrm{f}\left(\mathrm{\tau }\right)-\mathrm{f}\left(\mathrm{\tau }-\mathrm{\omega }\right)}{\mathrm{\omega }}\to {\mathrm{f}}^{\prime }\left(\mathrm{\tau }\right)$$E10

The above equation shows that the common derivative $\frac{d}{\mathit{d\tau }}$ is the limit of $D_{\omega }$ as $\omega \to 0$. That is,

$$D_{\omega }f\left(\tau \right)\le \left|\frac{\partial f}{\partial \tau }\right|=\underset{\omega \to 0}{lim}\frac{f\left(\tau +\omega \right)-f\left(\tau \right)}{\omega }$$E11

Now, the boundary element method (BEM) [47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80] is widely adopted for solving several engineering problems due to its easy implementation. In the BEM, only the boundary of the domain needs to be discretized, so it has a major advantage over other methods requiring full domain discretization [81, 82, 83, 84, 85, 86, 87] such as finite difference method (FDM), finite element method (FEM), and finite volume method (FVM) in engineering applications. 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. Previously scientists have proven that FEM covers more engineering applications than BEM which is more efficient for infinite domain problems. But currently BEM scientists have changed their thinking and vision on BEM, where the BEM researchers developed the BEM technique for solving inhomogeneous and nonlinear problems involving infinite and semi-infinite domains by using a lot of software like FastBEM and ExaFMM.

The main objective of this chapter is to introduce a novel memory-dependent derivative model for solving transient three-temperature nonlinear thermal stress problems in functionally graded anisotropic (FGA) smart structures. The governing equations of the considered model are nonlinear and very difficult if not impossible to solve analytically. Therefore, we develop a new efficient boundary element technique for solving such equations. Numerical results show the effects of MDD on the three-temperature distributions and the influence of MDD and anisotropy on the nonlinear thermal stresses of FGA smart structures. Also, numerical results demonstrate the validity and accuracy of the proposed model.

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 smart material problems, memory-dependent derivative history, and their applications. Section 2 describes the physical modeling of memory-dependent derivative problems of three-temperature nonlinear thermal stresses in FGA structures. Section 3 outlines the BEM implementation for obtaining the temperature field of the considered problem. Section 4 outlines the BEM implementation for obtaining the dispacement field of the considered problem. Section 5 introduces computing performance of the proposed model. Section 6 presents the new numerical results that describe the effects of memory-dependent derivative and anisotropy on the problem’s field variations. Lastly, Section 7 outlines the significant findings of this chapter.

2. Formulation of the problem

With reference to a Cartesian system $\left({\mathrm{x}}_1,{\mathrm{x}}_2,{\mathrm{x}}_3\right)$ with a configuration R bounded by a closed surface S as shown in Figure 1.

Figure 1.

Computational domain of the considered smart structure.

The governing equations for the transient three-temperature nonlinear thermal stresses problems of FGA smart structures with memory-dependent derivatives can be written as [35].

$${\sigma }_{\mathit{ij},j}+\rho F_i=\pi {\ddot{\mathrm{u}}}_i$$E12

$$D_{i,i}=0$$E13

where

$${\sigma }_{\mathit{ij}}={\left(x+1\right)}^m\left[C_{\mathit{\text{ijkl}}}e{\delta }_{\mathit{ij}}‐{\beta }_{\mathit{ab}}\left(T_{\alpha }‐T_{\alpha 0}+{\tau }_1{\dot{T}}_{\alpha }\right)\right]$$E14

$$D_i={\left(x+1\right)}^m\left[e_{\mathit{ijk}}{\epsilon }_{\mathit{jk}}+{\mathrm{f}}_{\mathit{ik}}{E}_k\right]{\epsilon }_{\mathit{ij}}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right),$$E15

where ${\sigma }_{\mathit{ij}}$, $F_i$, ${\epsilon }_{\mathit{ij}}$, ${\epsilon }_{\mathit{ijk}}$, $u_i$, and $\rho$ are the force stress tensor, mass force vector, strain tensor, alternate tensor, displacement vector, and density, respectively, $C_{\mathit{\text{ijkl}}}\left(C_{\mathit{\text{ijkl}}}=C_{\mathit{\text{klij}}}=C_{\mathit{\text{jikl}}}\right)$ is the constant elastic moduli, e is the dilatation, ${\beta }_{\mathit{ij}}\left({\beta }_{\mathit{ij}}={\beta }_{\mathit{ji}}\right)$ are the stress-temperature coefficients, $D_i$ is the electric displacement, $e_{\mathit{ijk}}$ is the piezoelectric tensor, ${\mathrm{f}}_{\mathrm{ik}}$ is the permittivity tensor, and $E_k$ is the electric field vector.

The two-dimensional three-temperature (2D-3T) radiative heat conduction equations can be expressed as

$$c_e\frac{\partial T_e\left(r,\tau \right)}{\partial \tau }-\frac{1}{\rho }\nabla \left[{\mathit{K}}_e\nabla T_e\left(r,\tau \right)\right]=-{\mathit{W}}_{\mathit{ei}}\left(T_e-T_i\right)-{\mathit{W}}_{\mathit{ep}}\left(T_e-T_p\right)$$E16

$$c_i\frac{\partial T_i\left(r,\tau \right)}{\partial \tau }-\frac{1}{\rho }\nabla \left[{\mathit{K}}_i\nabla T_i\left(r,\tau \right)\right]={\mathit{W}}_{\mathit{ei}}\left(T_e-T_i\right)$$E17

$$c_p{T}_p^3\frac{\partial T_p\left(r,\tau \right)}{\partial \tau }-\frac{1}{\rho }\nabla \left[{\mathit{K}}_p\nabla T_p\left(r,\tau \right)\right]={\mathit{W}}_{\mathit{ep}}\left(T_e-T_p\right)$$E18

where $e,i,\wedge p$ denote electron, ion, and phonon, respectively; $\left({\mathit{c}}_{\mathit{e}},{\mathit{c}}_{\mathit{i}},{\mathit{c}}_{\mathit{p}}\right)$, $\left({\mathit{K}}_{\mathit{e}},{\mathit{K}}_{\mathit{i}},{\mathit{K}}_{\mathit{p}}\right)$, and $\left({\mathit{T}}_{\mathit{e}},{\mathit{T}}_{\mathit{i}},{\mathit{T}}_{\mathit{p}}\right)$ are specific heat capacities, conductive coefficients, and temperature functions, respectively; ${\mathit{W}}_{\mathit{ei}}$ is the electron-ion coefficient; and ${\mathit{W}}_{\mathit{ep}}$ is the electron-phonon coefficient.

3. BEM solution of temperature field

This section concerns using a boundary element method to solve the temperature model.

The above 2D-3T radiative heat conduction Eqs. (16)-(18) can be expressed in the context of nonlinear thermal stresses of FGA smart structures as in [36].

$$\nabla \left[\left({\delta }_{1j}{\mathit{K}}_{\alpha }+{\delta }_{2j}{\mathit{K}}_{\alpha }^{\ast }\right)\nabla T_{\alpha }\left(r,\tau \right)\right]-\acute{\mathit{W}}\left(r,\tau \right)=c_{\alpha }{\mathit{\rho \delta }}_1{\delta }_{1j}{D}_{{\omega }_{\alpha }}{T}_{\alpha }\left(r,\tau \right)$$E19

which can be written in the following form:

$$L_{\mathit{ab}}{T}_{\alpha }\left(r,\tau \right)=f_{\mathit{ab}}$$E20

where

$$L_{\mathit{ab}}=\nabla \left[\left({\delta }_{1j}{\mathit{K}}_{\alpha }+{\delta }_{2j}{\mathit{K}}_{\alpha }^{\ast }\right)\nabla \right]$$E21

$$f_{\mathit{ab}}=\acute{\mathit{W}}\left(r,\tau \right)+\acute{\mathit{W}}\left(r,\tau \right)$$E22

where

$$\acute{\mathit{W}}\left(r,\tau \right)=\left\{\begin{array}{l}\rho {\mathit{W}}_{\mathit{ei}}\left(T_e-T_i\right)+\rho {\mathit{W}}_{\mathit{er}}\left(T_e-T_p\right)+\acute{\mathit{W}},\alpha =e,{\delta }_1=1\\ -\rho {\mathit{W}}_{\mathit{ei}}\left(T_e-T_i\right)+\acute{\mathit{W}},\alpha =i,{\delta }_1=1\\ -\rho {\mathit{W}}_{\mathit{er}}\left(T_e-T_p\right)+\acute{\mathit{W}},\alpha =p,{\delta }_1=T_p^3\end{array}\right.$$E23

$$\begin{array}{c}\acute{W}\left(r\text{,}\tau \right)=F\left(r\text{,}\tau \right)-\frac{{\delta }_{2j}{K}_{\alpha }}{{\omega }_{\alpha }}\underset{\tau -{\omega }_{\alpha }}{\overset{\tau }{{\displaystyle \int }}}K\left(\tau -\xi \right)\frac{\partial }{\partial \xi }\left({\nabla }^2{T}_{\alpha }\left(r\text{,}\tau \right)\right)d\xi \\ +\frac{\rho C_{\alpha }{\delta }_1{\delta }_{1j}}{{\omega }_{\alpha }}\underset{\tau -{\omega }_{\alpha }}{\overset{\tau }{{\displaystyle \int }}}K\left(\tau -\xi \right)\frac{\partial }{\partial \xi }\left(T_{\alpha }\left(r\text{,}\tau \right)\right)d\xi \\ +\frac{\rho C_{\alpha }\left({\tau }_0+{\delta }_{1j}{\tau }_2+{\delta }_{2j}\right)}{{\omega }_{\alpha }}\underset{\text{τ}-{\omega }_{\alpha }}{\overset{\text{τ}}{{\displaystyle \int }}}K\left(\tau -\xi \right)\frac{{\partial }^2}{\partial {\xi }^2}\left(T_{\alpha }\left(\text{r,}\tau \right)\right)d\xi \end{array}$$E24

$$F\left(r,\tau \right)={\beta }_{\mathit{ab}}{T}_{\alpha 0}\left[\stackrel{}{\AA }{\delta }_{1j}{\acute{u}}_{a,b}+\left({\tau }_0+{\delta }_{2j}\right){\acute{u}}_{a,b}\right]$$E25

and

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

where ${\delta }_{\mathit{ij}},\left(i,j=1,2\right)$, ${\omega }_{\alpha }\left(0\right)\left(\alpha =e,i\wedge p\right)$, and $K\left(\tau -\xi \right)$ are the Kronecker delta, delay times, and kernel function, respectively.

The total energy can be expressed as

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

Initial and boundary conditions can be expressed as

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

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

$${\mathbb{K}}_{\alpha }{\left.\frac{\partial T_{\alpha }}{\partial n}\right|}_{{\mathrm{\Gamma }}_2}=0,\alpha =e,i,p$$E30

By using the fundamental solutions $T_{\alpha }^{\ast }$ that satisfies the following differential equation:

$$L_{\mathit{ab}}{T}_{\alpha }^{\ast }=f_{\mathit{ab}}$$E31

Now, by implementing the technique of Fahmy [35], we can write (19) as

$${\mathit{CT}}_{\alpha }=\frac{D}{{\mathbb{K}}_{\alpha }}{\int }_O^{\tau }{\int }_S\left[T_{\alpha }{q}^{\ast }‐T_{\alpha }^{\ast }q\right]\mathit{dS}\mathit{d\tau }+\frac{D}{{\mathbb{K}}_{\alpha }}{\int }_O^{\tau }{\int }_{R}b{T}_{\alpha }^{\ast }\mathit{dR}\mathit{d\tau }+{\int }_R{\left.T_{\alpha }^i{T}_{\alpha }^{\ast }\right|}_{\tau =0}\mathit{dR}$$E32

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

$${\mathit{CT}}_{\alpha }={\int }_S\left[T_{\alpha }{q}^{\ast }-T_{\alpha }^{\ast }q\right]\mathit{dS}‐{\int }_R\frac{{\mathbb{K}}_{\alpha }}{D}\frac{\partial T_{\alpha }^{\ast }}{\partial \tau }{T}_{\alpha }\mathit{dR}$$E33

In order to transform the domain integral in (33) to the boundary, we approximate the temperature time derivative as

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

where $f^j\left(r\right)$ are known functions and $a^j\left(\tau \right)$ are unknown coefficients.

We assume that ${\stackrel{⌢}{T}}_{\alpha }^j$ is a solution of

$${\nabla }^2{\stackrel{⌢}{T}}_{\alpha }^j=f^j$$E35

Then, Eq. (33) leads to the following boundary integral equation

$$C{T}_{\alpha }={\int }_S\left[T_{\alpha }{q}^{\ast }‐T_{\alpha }^{\ast }q\right]\mathit{dS}+\sum _{i=1}^N{a}^j\left(\tau \right)D^{-1}\left(C{\stackrel{⌢}{T}}_{\alpha }^j‐{\int }_S\left[T_{\alpha }^j{q}^{\ast }‐{\widehat{q}}^j{T}_{\alpha }^{\ast }\right]\mathit{dS}\right)$$E36

where

$${\widehat{q}}^j=‐{\mathbb{K}}_{\alpha }\frac{\partial {\stackrel{⌢}{T}}_{\alpha }^j}{\partial n}$$E37

and

$$a^j\left(\tau \right)=\sum _{i=1}^N{f}_{\mathit{ji}}^{‐1}\frac{\partial T_{\alpha }\left(r_i,\tau \right)}{\partial \tau }$$E38

where $f_{\mathit{ji}}^{-1}$ are the coefficients of $F^{-1}$ which are defined as [58].

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

By discretizing Eq. (36) and using Eq. (38), we get [35].

$$C{\dot{T}}_{\alpha }+H{T}_{\alpha }=GQ$$E40

where Q is the heat flux vector and H and G are matrices.

The diffusion matrix can be defined as

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

where

$${\left\{\stackrel{⌢}{T}\right\}}_{\mathit{ij}}={\stackrel{⌢}{T}}^j\left(x_i\right)$$E42

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

To solve numerically Eq. (41), the functions $T_{\alpha }$ and q were interpolated as

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

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

where $0\le 0=\frac{\tau -{\tau }^m}{{\tau }^{m+1}‐{\tau }^m}\le 1$ determines the practical time $\tau$ of the current time step.

By time differentiation of Eq. (44), we obtain

$${\dot{T}}_{\alpha }=\frac{{\mathit{dT}}_{\alpha }}{\mathit{d\theta }}\frac{\mathit{d\theta }}{\mathit{d\tau }}=\frac{T_{\alpha }^{m+1}‐T_{\alpha }^m}{{\tau }^{m+1}‐{\tau }^m}=\frac{T_{\alpha }^{m+1}‐T_{\alpha }^m}{\mathrm{\Delta }{\tau }^m}$$E46

By substitution from (44)–(46) into (40), we get

$$\left(\frac{C}{\mathrm{\Delta }{\tau }^m}+\mathit{\theta H}\right)T_{\alpha }^{m+1}-\theta {\mathit{GQ}}^{m+1}=\left(\frac{C}{\mathrm{\Delta }{\tau }^m}‐\left(1-\theta \right)H\right)T_{\alpha }^m+\left(1-\theta \right){\mathit{GQ}}^m$$E47

By considering the initial and boundary conditions, we can write the following system of equations

$$\mathrm{a}X=\mathrm{b}$$E48

We apply an explicit staggered algorithm to solve the system (48) and obtain the temperature in terms of the displacement field.

4. BEM solution of displacement field

By using the weighted residual method, we can write (12) and (13) in the following form:

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

$${\int }_R\left(D_{,i}\right){\mathrm{\Phi }}_i^{\ast }\mathit{dR}=0$$E50

where

$$U_i=\rho F_i‐\rho {\ddot{u}}_i,$$E51

where $u_i^{\ast }$ and ${\mathrm{\Phi }}_i^{\ast }$ are weighting functions and $u_i$ and ${\mathrm{\Phi }}_i$ are approximate solutions.

Now, we assume the following boundary conditions:

$$u_i={\overline{u}}_i\mathit{on}{\mathrm{S}}_1$$E52

$${\lambda }_i={\sigma }_{\mathit{ij}}{n}_j={\overline{\lambda }}_i\mathit{on}{\mathrm{S}}_2$$E53

$$\mathrm{\Phi }=\overline{\mathrm{\Phi }}\mathit{on}{\mathrm{S}}_5$$E54

$$Q=\frac{\mathrm{\partial \Phi }}{\partial n}=\overline{Q}\mathit{on}{\mathrm{S}}_6$$E55

By integration by parts for the first term of Eqs. (49) and (50), we have

$$-{\int }_R{\sigma }_{\mathit{ij}}{u}_{i,j}^{\ast }\mathit{dR}+{\int }_R{U}_i{u}_i^{\ast }\mathit{dR}=-{\int }_{S_2}{\lambda }_i{u}_i^{\ast }\mathit{dS}$$E56

$$-{\int }_{R}D{\mathrm{\Phi }}_{i,i}^{\ast }\mathit{dR}=-{\int }_{S_6}{Q}_i{\mathrm{\Phi }}_i^{\ast }\mathit{dS}$$E57

Based on Huang and Liang [88], the boundary integral equation can be expressed as

$$\begin{array}{c}-\underset{R}{{\displaystyle \int }}{\sigma }_{ij,j}{u}_i^{\ast }dR+\underset{R}{{\displaystyle \int }}{U}_i{u}_i^{\ast }dR-\underset{R}{{\displaystyle \int }}D{\text{Φ}}_{i,i}^{\ast }dR=\underset{S_2}{{\displaystyle \int }}\left({\lambda }_i-{\overline{\lambda }}_i\right)u_i^{\ast }dS\\ +\underset{S_1}{{\displaystyle \int }}\left({\overline{u}}_i-u_i\right){\lambda }_i^{\ast }dS+\underset{S_6}{{\displaystyle \int }}\left(Q_i-{\overline{Q}}_i\right){\text{Φ}}_i^{\ast }dS\\ +\underset{S_5}{{\displaystyle \int }}\left({\overline{\text{Φ}}}_i-{\text{Φ}}_i\right)Q_i^{\ast }dS\end{array}$$E58

By integrating by parts for the left-hand side of (58), we get

$$\begin{array}{c}-\underset{R}{{\displaystyle \int }}{\sigma }_{ij}{\epsilon }_{ij}^{\ast }dR+\underset{R}{{\displaystyle \int }}{U}_i{u}_i^{\ast }dR-\underset{R}{{\displaystyle \int }}D{\text{Φ}}_{i,i}^{\ast }dR=-\underset{S_2}{{\displaystyle \int }}{\overline{\lambda }}_i{u}_i^{\ast }dS-\underset{S_1}{{\displaystyle \int }}{\lambda }_i{u}_i^{\ast }dS]\\ +\underset{S_1}{{\displaystyle \int }}\left({\overline{u}}_i-u_i\right){\lambda }_i^{\ast }dS-\underset{S_6}{{\displaystyle \int }}{\overline{Q}}_i{\text{Φ}}_i^{\ast }dS\\ -\underset{S_5}{{\displaystyle \int }}{Q}_i{\text{Φ}}_i^{\ast }dS+\underset{S_5}{{\displaystyle \int }}\left({\overline{\text{Φ}}}_i-{\text{Φ}}_i\right)Q_i^{\ast }dS\end{array}$$E59

Based on Eringen [89], the elastic stress can be expressed as

$${\sigma }_{\mathit{ij}}={\mathbb{A}}_{\mathit{\text{ijkl}}}{\epsilon }_{\mathit{kl}},$$E60

where

$${\mathbb{A}}_{\mathit{\text{ijkl}}}={\mathbb{A}}_{\mathit{\text{klij}}}$$E61

Hence, Eq. (59) can be rewritten as

$$\begin{array}{c}-\underset{R}{{\displaystyle \int }}{\sigma }^{ij\ast }{\epsilon }_{ij}dR+\underset{R}{{\displaystyle \int }}{U}_i{u}_i^{\ast }dR-\underset{R}{{\displaystyle \int }}D{\text{Φ}}_{i,i}^{\ast }dR=-\underset{S_2}{{\displaystyle \int }}{\overline{\lambda }}_i{u}_i^{\ast }dS-\underset{S_1}{{\displaystyle \int }}{\lambda }_i{u}_i^{\ast }dS\\ +\underset{S_1}{{\displaystyle \int }}\left({\overline{u}}_i-u_i\right){\lambda }_i^{\ast }dS-\underset{S_6}{{\displaystyle \int }}{\overline{Q}}_i{\text{Φ}}_i^{\ast }dS\\ -\underset{S_5}{{\displaystyle \int }}{Q}_i{\text{Φ}}_i^{\ast }dS+\underset{S_5}{{\displaystyle \int }}\left({\overline{\text{Φ}}}_i-{\text{Φ}}_i\right)Q_i^{\ast }dS\end{array}$$E62

By integration by parts again, we obtain

$${\int }_R{\sigma }_{\mathit{ij},i}^{\ast }{u}_i\mathit{dR}=-{\int }_S{u}_i^{\ast }{\lambda }_i\mathit{dS}-{\int }_S{\mathrm{\Phi }}_i^{\ast }{Q}_i\mathit{dS}+{\int }_S{\lambda }_i^{\ast }{u}_i\mathit{dS}+{\int }_S{Q}_i^{\ast }{\mathrm{\Phi }}_i\mathit{dS}$$E63

The weighting functions of $U_i={\mathrm{\Delta }}^n$ and $V_i=0$ along e₁ can be obtained as follows:

$${\sigma }_{1j,j}^{\ast }+{\mathrm{\Delta }}^n{e}_1=0$$E64

According to Dragos [90], the fundamental solution can be written as

$$u_i^{\ast }=u_{1i}^{\ast }{e}_1,{\mathrm{\Phi }}_i^{\ast }={\mathrm{\Phi }}_{1i}^{\ast }{e}_1,{\lambda }_i^{\ast }={\lambda }_{1i}^{\ast }{e}_1,{\mathrm{Q}}_i^{\ast }=Q_{1i}^{\ast }{e}_1$$E65

The weighting functions of $U_i=0$ and $V_i={\mathrm{\Delta }}^n$ along e₁ can be written as follows:

$${\sigma }_{\mathit{ij},j}^{\ast }=0$$E66

Based on Dragos [90], the fundamental solution can be obtained analytically as

$$u_i^{\ast }=u_{1i}^{\ast \ast }{e}_1,{\mathrm{\Phi }}_i^{\ast }={\mathrm{\Phi }}_{1i}^{\ast \ast }{e}_1,{\lambda }_i^{\ast }={\lambda }_{1i}^{\ast \ast }{e}_1,{\mathrm{Q}}_i^{\ast }=Q_{1i}^{\ast \ast }{e}_1$$E67

By using the weighting functions of (65) and (67) into (63), we have

$$C_{1i}^n{u}_i^n=-{\int }_S{\lambda }_{1i}^{\ast }{u}_i\mathit{dS}-{\int }_S{Q}_{1i}^{\ast }{\mathrm{\Phi }}_i\mathit{dS}+{\int }_S{u}_{1i}^{\ast }{\lambda }_i\mathit{dS}+{\int }_S{\mathrm{\Phi }}_{1i}^{\ast }{Q}_i\mathit{dS}$$E68

$$C_{1i}^n{\omega }_i^n=‐{\int }_S{\lambda }_{1i}^{\ast \ast }{u}_i\mathit{dS}-{\int }_S{Q}_{1i}^{\ast \ast }{\mathrm{\Phi }}_i\mathit{dS}+{\int }_S{u}_{1i}^{\ast \ast }{\lambda }_i\mathit{dS}+{\int }_S{\mathrm{\Phi }}_{1i}^{\ast \ast }{Q}_i\mathit{dS}$$E69

Thus, we can write

$$C^n{\mathrm{q}}^n=-{\int }_S{\mathrm{p}}^{\ast }\mathrm{q}\mathit{dS}+{\int }_S{\mathrm{q}}^{\ast }\mathrm{p}\mathit{dS}+{\int }_S{\mathrm{d}}^{\ast }\mathrm{\Phi }\mathit{ds}+{\int }_S{\mathrm{f}}^{\ast }\frac{\mathrm{\partial \Phi }}{\partial n}\mathit{dS}$$E70

where

$$\begin{array}{c}{C}^n=\left[\begin{array}{ll}{C}_{11}\hfill & C_{12}\hfill \\ C_{21}\hfill & C_{22}\hfill \end{array}\right],{\text{q}}^{\ast }=\left[\begin{array}{lll}{u}_{11}^{\ast }\hfill & u_{12}^{\ast }\hfill & 0\hfill \\ u_{21}^{\ast }\hfill & u_{22}^{\ast }\hfill & 0\hfill \\ u_{31}^{\ast \ast }\hfill & u_{32}^{\ast \ast }\hfill & 0\hfill \end{array}\right],{\text{p}}^{\ast }=\left[\begin{array}{lll}{\lambda }_{11}^{\ast }\hfill & {\lambda }_{12}^{\ast }\hfill & 0\hfill \\ {\lambda }_{21}^{\ast }\hfill & {\lambda }_{22}^{\ast }\hfill & 0\hfill \\ {\lambda }_{31}^{\ast \ast }\hfill & {\lambda }_{32}^{\ast \ast }\hfill & 0\hfill \end{array}\right],\text{q}=\left[\begin{array}{l}{u}_1\hfill \\ u_2\hfill \\ {\omega }_3\hfill \end{array}\right],\\ \text{p}=\left[\begin{array}{l}{\lambda }_1\hfill \\ \lambda \hfill \\ {\mu }_3\hfill \end{array}\right],{\text{d}}^{\ast }=\left[\begin{array}{c}{\text{d}}_1^{\ast }\\ {\text{d}}_2^{\ast }\\ 0\end{array}\right],{\text{f}}^{\ast }=\left[\begin{array}{c}{\text{f}}_1^{\ast }\\ {\text{f}}_2^{\ast }\\ 0\end{array}\right]\end{array}$$E71

In order to solve (70) numerically, we suppose the following definitions:

$$\mathrm{q}=\psi {\mathrm{q}}^j,\mathrm{p}=\psi {\mathrm{p}}^j,\mathrm{\Phi }={\psi }_0{\mathrm{\Phi }}^j,\frac{\mathrm{\partial \Phi }}{\partial n}={\psi }_0{\left(\frac{\mathrm{\partial \Phi }}{\partial n}\right)}^j$$E72

Substituting from (72) into (70) and discretizing the boundary, we obtain

$$C^n{\mathrm{q}}^n=\sum _{j=1}^{N_e}\left[-{\int }_{{\mathrm{\Gamma }}_j}{\mathrm{p}}^{\ast }\psi d\mathrm{\Gamma }\right]{\mathrm{q}}^j+\sum _{j=1}^{N_e}\left[{\int }_{{\mathrm{\Gamma }}_j}{\mathrm{q}}^{\ast }\psi d\mathrm{\Gamma }\right]{\mathrm{p}}^j+\sum _{j=1}^{N_e}\left[{\int }_{{\mathrm{\Gamma }}_j}{\mathrm{d}}^{\ast }{\psi }_{0}d\mathrm{\Gamma }\right]{\mathrm{\Phi }}^j+\sum _{j=1}^{N_e}\left[{\int }_{{\mathrm{\Gamma }}_j}{\mathrm{f}}^{\ast }{\psi }_{0}d\mathrm{\Gamma }\right]{\left(\frac{\mathrm{\partial \Phi }}{\partial n}\right)}^j$$E73

Equation after integration can be written as

$$C^i{\mathrm{q}}^i=-\sum _{j=1}^{N_e}{\widehat{\mathbb{H}}}^{\mathit{ij}}{\mathrm{q}}^j+\sum _{j=1}^{N_e}{\widehat{\mathbb{G}}}^{\mathit{ij}}{\mathrm{p}}^j+\sum _{j=1}^{N_e}{\widehat{\mathbb{D}}}^{\mathit{ij}}{\mathrm{\Phi }}^i+\sum _{j=1}^{N_e}{\widehat{\mathbb{F}}}^{\mathit{ij}}{\left(\frac{\mathrm{\partial \Phi }}{\partial n}\right)}^j$$E74

By using the following representation:

$${\mathbb{H}}^{\mathit{ij}}=\left\{\begin{array}{ll}{\widehat{\mathbb{H}}}^{\mathit{ij}}& \mathit{\text{if}}i\ne j\\ {\widehat{\mathbb{H}}}^{\mathit{ij}}+C^i& \mathit{\text{if}}i=j\end{array}\right.$$E75

Thus, we can write (74) as follows:

$$\sum _{j=1}^{N_e}{\mathbb{H}}^{\mathit{ij}}{\mathrm{q}}^j=\sum _{j=1}^{N_e}{\widehat{\mathbb{G}}}^{\mathit{ij}}{\mathrm{p}}^j+\sum _{j=1}^{N_e}{\widehat{\mathbb{D}}}^{\mathit{ij}}{\mathrm{\Phi }}^j+\sum _{j=1}^{N_e}{\widehat{\mathbb{F}}}^{\mathit{ij}}{\left(\frac{\mathrm{\partial \Phi }}{\partial n}\right)}^j$$E76

The global matrix equation for all $\mathbf{i}$ nodes can be expressed as

$$\mathbb{HQ}=\mathbb{G}\mathbb{P}+\mathbb{D}\mathrm{\Theta }+\mathbb{FS}$$E77

where $\mathbb{Q}$ is the displacement vector, $\mathbb{P}$ is the traction vector, $\mathrm{\Theta }$ is the electric potential vector, and $\mathbb{S}$ is the electric potential gradient vector.

Substituting the boundary conditions into (77), we obtain the following system of equations:

$$\mathbb{AX}=\mathbb{B}$$E78

We apply an explicit staggered algorithm to solve the system (78) and obtain the temperature and displacement fields as follows:

1. From Eq. (48) we obtain the temperature field in terms of the displacement field.

2. We predict the displacement field and solve the resulted equation for the temperature field.

3. We correct the displacement field using the computed temperature field for Eq. (78).

An explicit staggered algorithm based on communication-avoiding Arnoldi as described in Hoemmen [91] is very suitable for efficient implementation in Matlab (R2019a) with the aim of specifically improving its performance for the solution of the resulting linear algebraic systems.

5. Computational performance of the problem

According to Fahmy [35], the computer performance with simulation can be computed based on account and communication process, elements underlying the hardware and functional computation. The main objective of our proposed technique during simulation process is to use the preconditioners which are efficient to improve the overall CPU utilization of the cluster, accelerate the iterative method, and reduce the input/output and the interprocessor communication costs. Also, Fahmy [35] compared the communication-avoiding Krylov methods that are based on the s-step Krylov methods such as communication-avoiding generalized minimal residual (CA-GMRES) of Saad and Schultz [92], communication-avoiding Arnoldi (CA-Arnoldi) of the Arnoldi [93] and communication-avoiding Lanczos (CA-Lanczos) of Lanczos [94], with their corresponding standard Krylov methods. CA-Arnoldi which is also called Arnoldi (s, t) algorithm is different from standard Arnoldi (s) $\left(s,t=1\right)$, where s is the number of inner iteration steps and t is the number of outer iteration steps. According to [35], the CA-Arnoldi has numerical stability, convergence, and performance due to the implementation of algorithm shown in Figure 2, which is based on the QR factorization update and block classical Gram-Schmidt (block CGS) approach or block modified Gram-Schmidt (block MGS) approach where

Figure 2.

CA-Arnoldi iteration algorithm.

$$V_k=\left[v_{\mathit{sk}+1},v_{\mathit{sk}+2},\dots ,v_{\mathit{sk}+s}\right]$$E79

and

$${\mathcal{Q}}_k=\left[{\mathrm{Q}}_0,{\mathrm{Q}}_1,\dots ,{\mathrm{Q}}_{k-1}\right]$$E80

The generalized minimal residual (GMRES) method of Saad and Schultz [92] is a Krylov subspace method for solving nonsymmetric linear systems. The CA-GMRES algorithm is based on Arnoldi (s, t) and equivalent to standard GMRES in exact arithmetic. Also, the GMRES or CA-GMRES are convergent at the same rate for problems, but Hoemmen [91] proved that CA-GMRES algorithm shown in Figure 3 converges for the s-step basis lengths and restart lengths used for obtaining maximum performance. Lanczos method can be considered as a special case of Arnoldi method for symmetric and real case of A or Hermitian and complex case of A. Symmetric Lanczos which is also called Lanczos is different from nonsymmetric Lanczos. We implemented a communication-avoiding version of symmetric Lanczos (CA-Lanczos) for solving symmetric positive definite (SPD) eigenvalue problems. Also, we implement CA-Lanczos iteration algorithm shown in Figure 4, which is also called Lanczos (s, t), where s is the s-step basis length and t is the outer iterations number before restart. This algorithm is based on using rank revealing-tall skinny QR-block Gram-Schmidt (RR-TSQR-BGS) orthogonalization method which connects between TSQR and block Gram-Schmidt, where we have been using the right-shifted basis matrix at outer iteration k as follows:

Figure 3.

CA-GMRES iteration algorithm.

Figure 4.

CA-Lanczos iteration algorithm.

$$V_k^{\prime }=\left[V_{\mathit{sk}+2},\dots ,v_{\mathit{sk}+s}\right]$$E81

and

$$V_k^{\prime }=\left[V_k^{\prime },v_{\mathit{sk}+s+1}\right]$$E82

For more details about the considered preconditioners and algorithms, we refer the interested readers to [91].

The main objective of this section is to implement an accurate and robust preconditioning technique for solving the dense nonsymmetric algebraic system of linear equations arising from the BEM. So, a communication-avoiding Arnoldi of the Arnoldi [93] has been implemented for solving the resulting linear systems in order to reduce the iteration number and CPU time. The BEM discretization is employed in 1280 quadrilateral elements, with 3964 degrees of freedom (DOF). A comparative performance of preconditioned Krylov subspace solvers (CA-Arnoldi, CA-GMRES, and CA-Lanczos) has been shown in Table 1, where the number of DOF is 3964 and “–” was defined as the divergence process. From the results of Table 1. The CA-Arnoldi, CA-GMRES, and CA-Lanczos are more cost-effective than the other Krylov subspace methods Arnoldi, GMRES, and Lanczos, respectively. Also, CA-Arnoldi, CA-GMRES, and CA-Lanczos have been compared with each other in Table 2. It can be seen from this table that the performance of CA-Arnoldi is superior than the other iterative methods.

Methods	Preconditioning techniques	Iterations	Residual	Time of each iterative step (s)	Time of solution
Direct methods	NO	—	—	—	9 min 50 s
Arnoldi	NO	174	7.21E–07	3.85	11 min 25 s
	JOBI	26	5.22E–07	3.86	2 min 38 s
	BJOB	22	1.34E–06	3.86	2 min 23 s
	ILU3	47	1.66E–06	3.84	4 min 2 s
	ILU5	48	1.38E–06	3.89	4 min 6 s
	DILU	48	1.53E–06	5.45	4 min 18 s
CA–Arnoldi	NO	360	6.96E–07	1.95	11 min 53 s
	JOBI	20	4.42E–07	1.96	1 min 30 s
	BJOB	20	2.30E–08	1.96	1 min 30 s
	ILU3	40	7.87E–07	1.96	2 min 11 s
	ILU5	60	1.28E–08	1.96	2 min 48 s
	DILU	60	1.59E–07	3.07	4 min 1 s
GMRES	NO	280	2.36E–08	1.90	6 min 20 s
	JOBI	40	5.01E–13	1.91	2 min 10 s
	BJOB	40	2.05E–11	1.91	2 min 10 s
	ILU3	40	4.70E–08	1.91	2 min 10 s
	ILU5	40	3.13E–08	2.60	2 min 10 s
	DILU	40	6.19E–08	3.07	2 min 48 s
CA–GMRES	NO	120	6.89E–07	3.78	7 min 57 s
	JOBI	12	1.00E–05	3.76	1 min 41 s
	BJOB	12	2.22E–06	3.76	1 min 42 s
	ILU3	26	3.63E–06	3.75	2 min 34 s
	ILU5	22	4.05E–06	3.75	2 min 20 s
	DILU	25	5.19E–06	5.93	3 min 18 s
Lanczos	NO	135	7.24E–07	3.80	8 min 41 s
	JOBI	22	4.87E–07	3.75	2 min 33 s
	BJOB	18	9.27E–07	5.18	3 min 2 s
	ILU3	42	2.41E–07	3.81	3 min 48 s
	ILU5	36	6.41E–07	3.78	3 min 18 s
	DILU	38	2.04E–07	5.00	3 min 32 s
CA–Lanczos	NO	129	1.30E–04	3.75	9 min 22 s
	JOBI	16	8.64E–07	3.76	2 min 3s
	BJOB	14	1.69E–07	3.77	2 min 0 s
	ILU3	24	9.29E–07	3.87	2 min 31 s
	ILU5	31	1.91E–07	3.90	3 min 1 s
	DILU	27	8.11E–07	5.95	3 min 31 s

Table 1.

Performances of preconditioned Krylov subspace iterative methods for DOF 3964.

Solvers		DOF
		965	1505	3380	3964	6005
CA–Arnoldi	Residual	6.81E–12	5.38E–12	4.13E–11	4.17E–11	7.57E–11
	CPU time (s)	4.96	10.78	99.24	134.26	293.29
	Iterations	25	25	25	25	25
CA–GMRES	Residual	2.98E–12	1.90E–12	1.28E–11	1.36E–11	1.22E–11
	CPU time (s)	5.06	11.49	126.38	164.09	445.51
	Iterations	50	50	50	50	50
CA– Lanczos	Residual	7.20E–11	3.35E–11	2.72E–11	3.97E–11	8.33E–11
	CPU time (s)	5.05	11.47	139.07	180.49	514.72
	Iterations	22	26	28	30	32

Table 2.

The CPU time and the number of iterations for some communication–avoiding Krylov subspace solvers.

6. Numerical results and discussion

In order to illustrate the numerical results of this study, we consider a monoclinic graphite-epoxy as an anisotropic smart material which has the following constants [35].

The elasticity tensor is expressed as

$$C_{p\mathit{jkl}}=\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}$$E83

The mechanical temperature coefficient is

$${\beta }_{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$$E84

The thermal conductivity tensor is

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

Mass density $\rho =7820\mathrm{kg}/{\mathrm{m}}^2$ and heat capacity $c=461\mathrm{J}/\mathrm{kg}\mathrm{k}$.

The technique that has been proposed in the current chapter can be applicable to a wide range of three-temperature nonlinear thermal stress problems of FGA structures. The main aim of this chapter is to assess the impact of MDD and anisotropy on the three-temperature nonlinear thermal stress distributions.

The proposed technique that has been implemented in the current study can be applicable to a wide variety of FGA smart structure problems involving three temperatures. All the physical parameters satisfy the initial and boundary conditions. The efficiency of our BEM modeling technique has been improved using an explicit staggered algorithm based on communication-avoiding Arnoldi procedure to decrease the computation time.

Figure 5 shows the variations of the three temperatures ${\mathrm{T}}_{\mathrm{e}},{\mathrm{T}}_{\mathrm{i}}\text{and}{\mathrm{T}}_{\mathrm{p}}$ with the time $\mathrm{\tau }$ in the presence of MDD. Figure 6 shows the variations of the three temperatures ${\mathrm{T}}_{\mathrm{e}},{\mathrm{T}}_{\mathrm{i}}\text{and}{\mathrm{T}}_{\mathrm{p}}$ with the time $\mathrm{\tau }$ in the presence of MDD. It can be seen from Figures 5 and 6 that the MDD has a significant effect on the temperature distributions.

Figure 5.

Variation of the three-temperature (with memory) with time τ.

Figure 6.

Variation of the three-temperature (without memory) with time τ.

In order to study the anisotropy and MDD effects on the nonlinear thermal stresses, we assume the following four cases: A, B, C, and D, where case A denotes the nonlinear thermal stress distribution in the isotropic material without MDD effect, case B denotes the nonlinear thermal stress distribution in isotropic material with MDD effect, case C denotes the nonlinear thermal stress distribution in anisotropic material without MDD effect, and case D denotes nonlinear thermal stress distribution in anisotropic material with MDD effect.

Figures 7–9 show the variation of the nonlinear thermal stresses ${\mathrm{\sigma }}_{11},{\mathrm{\sigma }}_{12}\text{and}{\mathrm{\sigma }}_{22}$ with the time $\mathrm{\tau }$. It is clear from these figures that both anisotropy and MDD have a significant influence on the nonlinear thermal stress distributions.

Figure 7.

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

Figure 8.

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

Figure 9.

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

Since there are no available results for the considered problem in the literature. Therefore, we only considered the one-dimensional special case for the variations of the nonlinear thermal stress σ₁₁ with the time $\mathrm{\tau }$ as shown in Figure 10. The validity and accuracy of our proposed technique was confirmed by comparing graphically our BEM results with those obtained using the FDM of Pazera and Jędrysiak [95] and FEM of Xiong and Tian [96] results based on replacing one-temperature heat conduction with the total three-temperature T $\left(\mathrm{T}={\mathrm{T}}_{\mathrm{e}}+{\mathrm{T}}_{\mathrm{i}}+{\mathrm{T}}_{\mathrm{r}}\right)$ heat conduction. It can be noticed that the BEM results are found to agree very well with the FDM and FEM results.

Figure 10.

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

7. Conclusion

The main aim of this chapter is to introduce a new MDD model based on BEM for obtaining the transient three-temperature nonlinear thermal stresses in FGA smart structures. The governing equations of this model are very hard to solve analytically because of nonlinearity and anisotropy. To overcome this, we propose a new boundary element formulation for solving such equations. Since the CA kernels of the s-step Krylov methods are faster than the kernels of standard Krylov methods. Therefore, we used an explicit staggered algorithm based on CA-Arnoldi procedure to solve the resulted linear equations. The computational performance of the proposed technique has been performed using communication-avoiding Arnoldi procedure. The numerical results are presented highlighting the effects of MDD on the three-temperature distributions and the influence of MDD and anisotropy on the nonlinear thermal stresses of FGA smart structures. The numerical results also demonstrate the validity and accuracy of the proposed technique. It can be concluded from numerical results of our current general problem that all generalized and nonlinear generalized thermoelasticity theories can be combined with the three-temperature radiative heat conduction to describe the deformation of FGA smart structures in the context of memory-dependent derivatives. From the research that has been performed, it is possible to conclude that the proposed BEM technique is effective and stable for transient three-temperature thermal stress problems in FGA smart structures.

The numerical results for our complex and general problem can provide data references for computer scientists and engineers, geotechnical and geothermal engineers, designers of new materials, and researchers in material science as well as for those working on the development of anisotropic smart structures. In the application of three-temperature theories in advanced manufacturing technologies, with the development of soft machines and robotics in biomedical engineering and advanced manufacturing, transient thermal stresses 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 smart materials and structures.