Boundary Element Mathematical Modelling and Boundary Element Numerical Techniques for Optimization of Micropolar Thermoviscoelastic Problems in Solid Deformable Bodies

Abstract

The main objective of this chapter is to introduce a new theory called three-temperature nonlinear generalized micropolar thermoviscoelasticity. Because of strong nonlinearity of simulation and optimization problems associated with this theory, the numerical solutions for problems related with the proposed theory are always very difficult and require the development of new numerical techniques. So, we propose a new boundary element technique for simulation and optimization of such problems based on genetic algorithm (GA), free form deformation (FFD) method and nonuniform rational B-spline curve (NURBS) as the shape optimization technique. In the formulation of the considered problem, the profiles of the considered objects are determined by FFD method, where the FFD control points positions are treated as genes, and then the chromosomes profiles are defined with the genes sequence. The population is founded by a number of individuals (chromosomes), where the objective functions of individuals are determined by the boundary element method (BEM). The numerical results verify the validity and accuracy of our proposed technique.

1. Introduction

The classical thermo-elasticity (CTE) theory which was introduced by Duhamel [1] and Neumann [2] characterized the strain-temperature gradients equations in an elastic body, but it has two shortcomings contrary to physical observations: First, the heat conduction equation of this theory does not include any elastic terms. Second, the heat conduction equation is of a parabolic type predicting infinite speeds of thermal waves. Biot [3] developed the classical coupled thermo-elasticity (CCTE) theory to overcome the first shortcoming in CTE. However, both theories share the second shortcoming. So, several generalized thermoelasticity theories that predict finite speeds of propagation for heat waves have been developed such as extended thermo-elasticity (ETE) theory of Lord and Shulman [4], temperature-rate-dependent thermo-elasticity (TRDTE) theory of Green and Lindsay [5], three linear generalized thermoelasticity theories (type I, II and III) of Green and Naghdi [6, 7], low-temperature thermoelasticity (LTTE) model of Hetnarski and Ignaczak [8], the dual phase-lag (DPL) heat conduction equation of Tzou [9, 10] which has been developed taking into consideration the phonons-electrons interactions to obtain dual phase-lag thermoelasticity (DPLTE) [11, 12], and three-phase-lag thermoelasticity (TPLTE) model of Choudhuri [13] who takes into consideration the phase-lags of heat flux, temperature gradient and thermal displacement gradient. Chen and Gurtin [14], introduced the theory of two-temperatures (conductive temperature and thermodynamic temperature) heat conduction in the context of elastic bodies, then Youssef [15] extended this theory to generalized thermoelasticity. Fahmy [16] introduced three-temperature nonlinear generalized micropolar-magneto-thermoelasticity theory and developed a new boundary element technique for Modeling and Simulation of complex problems associated with this theory. Theory of micropolar elasticity [17, 18] has been developed for studying the mechanical behavior of polymers and elastomers and applied in many applications [19, 20, 21, 22, 23, 24]. Then, Eringen [25] and Nowacki [26] extended it to micropolar thermoelasticity, and then implemented in various applications [27, 28, 29]. Because of strong nonlinearity of three-temperatures radiative heat conduction equations, the numerical solution and simulation of such problems are always difficult and require the development of new numerical schemes [30, 31]. In comparison with other numerical methods [32, 33, 34], the boundary element method has been successfully applied and performed for solving various applications [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60]. The boundary element technique has been formulated in the context of micropolar thermoelasticity by Sladek and Sladek [61, 62, 63] and Huang and Liang [64]. Through the current paper, the term three-temperatures introduced for the first time in the field of nonlinear generalized micropolar thermoviscoelasticity. Recently, evolutionary algorithms [65, 66] have received much attention of researchers. The genetic algorithm (GA) can deal with the multi-objective and complex shapes problems. Also, it could reach an optimal solution with highly reduced computational cost.

The main aim of this article is to introduce a new theory called nonlinear generalized micropolar thermoviscoelasticity involving three temperatures. Because of strong nonlinearity, it is very difficult to solve the problems related to this theory analytically. Therefore, we propose a new boundary element model for simulation and optimization of three temperatures nonlinear generalized micropolar thermoviscoelastic problems associated with this theory. The genetic algorithm (GA) was implemented based on free form deformation (FFD) technique and nonuniform rational B-spline (NURBS) curve as an optimization technique for the considered problems. The numerical results demonstrate the validity and accuracy of our proposed model.

2. Formulation of the problem

The governing equations for three-temperature anisotropic generalized micropolar thermoviscoelasticity problems can be expressed as [58]

$${\mathrm{\sigma }}_{\mathrm{ij},\mathrm{j}}+\mathrm{\rho }{\mathrm{F}}_{\mathrm{i}}=\mathrm{\rho }{\ddot{\mathrm{u}}}_{\mathrm{i}}$$ E1

$${\mathrm{m}}_{\mathrm{ij},\mathrm{j}}+{\mathrm{\epsilon }}_{\mathrm{ijk}}{\mathrm{\sigma }}_{\mathrm{jk}}+\mathrm{\rho }{\mathrm{M}}_{\mathrm{i}}=\mathrm{J\rho }{\ddot{\mathrm{\omega }}}_{\mathrm{i}}$$ E2

where

$${\mathrm{\sigma }}_{\mathrm{ij}}={\mathrm{C}}_{\text{ijkl}}\mathrm{\aleph }\mathrm{e}{\mathrm{\delta }}_{\mathrm{ij}}+\mathrm{\alpha }ˇ\left({\mathrm{u}}_{\mathrm{j},\mathrm{i}}-{\mathrm{\epsilon }}_{\mathrm{ijk}}{\mathrm{\omega }}_{\mathrm{k}}\right)-{\mathrm{\beta }}_{\mathrm{ij}}{\mathrm{T}}_{\mathrm{\alpha }}\left({\mathrm{C}}_{\text{ijkl}}={\mathrm{C}}_{\text{klij}}={\mathrm{C}}_{\text{jikl}},{\mathrm{\beta }}_{\mathrm{ij}}={\mathrm{\beta }}_{\mathrm{ji}}\right)$$ E3

$${\mathrm{m}}_{\mathrm{ij}}=\mathrm{\alpha }{\mathrm{\omega }}_{\mathrm{k},\mathrm{k}}{\mathrm{\delta }}_{\mathrm{ij}}+\overline{\mathrm{\alpha }}{\mathrm{\omega }}_{\mathrm{i},\mathrm{j}}+\mathrm{\alpha }̿{\mathrm{\omega }}_{\mathrm{j},\mathrm{i}}$$ E4

$${\mathrm{ϵ}}_{\mathrm{ij}}={\mathrm{\epsilon }}_{\mathrm{ij}}-{\mathrm{\epsilon }}_{\mathrm{ijk}}\left({\mathrm{r}}_{\mathrm{k}}-{\mathrm{\omega }}_{\mathrm{k}}\right),{\mathrm{\epsilon }}_{\mathrm{ij}}=\frac{1}{2}\left({\mathrm{u}}_{\mathrm{i},\mathrm{j}}+{\mathrm{u}}_{\mathrm{j},\mathrm{i}}\right),{\mathrm{r}}_{\mathrm{i}}=\frac{1}{2}{\mathrm{\epsilon }}_{\mathrm{ikl}}{\mathrm{u}}_{\mathrm{l},\mathrm{k}}$$ E5

The two-dimension three-temperature (2D-3T) radiative heat conduction equations can be expressed as [53]

$${\mathrm{c}}_{\mathrm{e}}\frac{\partial {\mathrm{T}}_{\mathrm{e}}\left(\mathrm{r},\mathrm{\tau }\right)}{\mathrm{\partial \tau }}-\frac{1}{\mathrm{\rho }}\nabla \left[{\mathbb{K}}_{\mathrm{e}}\nabla {\mathrm{T}}_{\mathrm{e}}\left(\mathrm{r},\mathrm{\tau }\right)\right]=-{\mathbb{W}}_{\mathrm{ei}}\left({\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{i}}\right)-{\mathbb{W}}_{\mathrm{ep}}\left({\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{p}}\right)$$ E6

$${\mathrm{c}}_{\mathrm{i}}\frac{\partial {\mathrm{T}}_{\mathrm{i}}\left(\mathrm{r},\mathrm{\tau }\right)}{\mathrm{\partial \tau }}-\frac{1}{\mathrm{\rho }}\nabla \left[{\mathbb{K}}_{\mathrm{i}}\nabla {\mathrm{T}}_{\mathrm{i}}\left(\mathrm{r},\mathrm{\tau }\right)\right]={\mathbb{W}}_{\mathrm{ei}}\left({\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{i}}\right)$$ E7

$$\frac{4}{\mathrm{\rho }}{\mathrm{c}}_{\mathrm{p}}{\mathrm{T}}_{\mathrm{p}}^3\frac{\partial {\mathrm{T}}_{\mathrm{p}}\left(\mathrm{r},\mathrm{\tau }\right)}{\mathrm{\partial \tau }}-\frac{1}{\mathrm{\rho }}\nabla \left[{\mathbb{K}}_{\mathrm{p}}\nabla {\mathrm{T}}_{\mathrm{p}}\left(\mathrm{r},\mathrm{\tau }\right)\right]={\mathbb{W}}_{\mathrm{ep}}\left({\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{p}}\right)$$ E8

3. A new mathematical modelling of nonlinear generalized micropolar thermoviscoelasticity problem

With reference to a Cartesian coordinate system $\left(x_1,x_2,x_3\right)$ , we consider an anisotropic micropolar thermoviscoelastic structure occupies the region $\mathrm{R}$ which bounded by a closed surface $\mathrm{S}$ , and ${\mathrm{S}}_{\mathrm{i}}\left(\mathrm{i}=1,2,3,4\right)$ denotes subsets of $\mathrm{S}$ such that ${\mathrm{S}}_1+{\mathrm{S}}_2={\mathrm{S}}_3+{\mathrm{S}}_4=\mathrm{S}.$

3.1 BEM simulation for temperature field

The 2D-3T radiative heat conduction Eqs. (6)–(8) can be expressed as [53]

$$\nabla \left[{\mathbb{K}}_{\mathrm{\alpha }}\nabla {\mathrm{T}}_{\mathrm{\alpha }}\left(\mathrm{r},\mathrm{\tau }\right)\right]+\overline{\mathbb{W}}\left(\mathrm{r},\mathrm{\tau }\right)={\mathrm{c}}_{\mathrm{\alpha }}\mathrm{\rho }{\mathrm{\delta }}_1\frac{\partial {\mathrm{T}}_{\mathrm{\alpha }}\left(\mathrm{r},\mathrm{\tau }\right)}{\mathrm{\partial \tau }}$$ E9

where

$$\overline{\mathbb{W}}\left(\mathrm{r},\mathrm{\tau }\right)=\left\{\begin{array}{l}-\mathrm{\rho }{\mathbb{W}}_{\mathrm{ei}}\left({\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{i}}\right)–\mathrm{\rho }{\mathbb{W}}_{\mathrm{ep}}\left({\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{p}}\right),\mathrm{\alpha }=\mathrm{e},{\mathrm{\delta }}_1=1\\ \mathrm{\rho }{\mathbb{W}}_{\mathrm{ei}}\left({\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{i}}\right),\mathrm{\alpha }=\mathrm{i},{\mathrm{\delta }}_1=1\\ \mathrm{\rho }{\mathbb{W}}_{\mathrm{ep}}\left({\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{p}}\right),\mathrm{\alpha }=\mathrm{p},{\mathrm{\delta }}_1=\frac{4}{\mathrm{\rho }}{\mathrm{T}}_{\mathrm{p}}^3\end{array}\right.$$ E10

and

$${\mathbb{W}}_{\mathrm{ei}}=\mathrm{\rho }{\mathbb{A}}_{\mathrm{ei}}{\mathrm{T}}_{\mathrm{e}}^{-2/3},{\mathbb{W}}_{\mathrm{ep}}=\mathrm{\rho }{\mathbb{A}}_{\mathrm{ep}}{\mathrm{T}}_{\mathrm{e}}^{-1/2},{\mathbb{K}}_{\mathrm{\alpha }}={\mathbb{A}}_{\mathrm{\alpha }}{\mathrm{T}}_{\mathrm{\alpha }}^{5/2},\mathrm{\alpha }=\mathrm{e},\mathrm{i},{\mathbb{K}}_{\mathrm{p}}={\mathbb{A}}_{\mathrm{p}}{\mathrm{T}}_{\mathrm{p}}^{3+\mathbb{B}}$$ E11

where parameters ${\mathrm{c}}_{\mathrm{\alpha }}$ , ${\mathbb{A}}_{\mathrm{\alpha }}\left(\mathrm{\alpha }=\mathrm{e},\mathrm{i},\mathrm{p}\right)$ , $\mathbb{B}$ , ${\mathbb{A}}_{\mathrm{ei}}$ , ${\mathbb{A}}_{\mathrm{ep}}$ are constant inside each subdomain, but they are discontinuous on the interfaces between subdomains.

The total energy of unit mass can be described by

$$\mathrm{P}={\mathrm{P}}_{\mathrm{e}}+{\mathrm{P}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{p}},{\mathrm{P}}_{\mathrm{e}}={\mathrm{c}}_{\mathrm{e}}{\mathrm{T}}_{\mathrm{e}},{\mathrm{P}}_{\mathrm{i}}={\mathrm{c}}_{\mathrm{i}}{\mathrm{T}}_{\mathrm{i}},{\mathrm{P}}_{\mathrm{p}}=\frac{1}{\mathrm{\rho }}{\mathrm{c}}_{\mathrm{p}}{\mathrm{T}}_{\mathrm{p}}^4$$ E12

Initial and boundary conditions can be written as

$${\mathrm{T}}_{\mathrm{\alpha }}\left(\mathrm{x},\mathrm{y},0\right)={\mathrm{T}}_{\mathrm{\alpha }}^0\left(\mathrm{x},\mathrm{y}\right)={\mathrm{g}}_1\left(\mathrm{x},\mathrm{\tau }\right)$$ E13

$${\left.{\mathbb{K}}_{\mathrm{\alpha }}\frac{\partial {\mathrm{T}}_{\mathrm{\alpha }}}{\mathrm{\partial n}}\right|}_{{\mathrm{\Gamma }}_1}=0,\mathrm{\alpha }=\mathrm{e},\mathrm{i},{\left.{\mathrm{T}}_{\mathrm{p}}\right|}_{{\mathrm{\Gamma }}_1}={\mathrm{g}}_2\left(\mathrm{x},\mathrm{\tau }\right)$$ E14

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

we use the time-dependent fundamental solution which is a solution of the following differential equation

$$\mathrm{D}{\nabla }^2{\mathrm{T}}_{\mathrm{\alpha }}+\frac{\partial {\mathrm{T}}_{\mathrm{\alpha }}^{\ast }}{\mathrm{\partial n}}=-\mathrm{\delta }\left(\mathrm{r}-{\mathrm{p}}_{\mathrm{i}}\right)\mathrm{\delta }\left(\mathrm{\tau }-\mathrm{r}\right),\mathrm{D}=\frac{{\mathbb{K}}_{\mathrm{\alpha }}}{\mathrm{\rho c}}$$ E16

In which the points ${\mathrm{p}}_{\mathrm{i}}$ are the singularities, where the temperatures are not defined there. Singular integrals are those whose kernels are not defined at the singularities on the integration domain $\mathrm{R}$ . They are defined by eliminating a small space including the singularity, and obtaining the limit when this small space tends to zero [40, 46].

The boundary integral equation corresponding to our considered heat conduction can be written as in Fahmy [46, 47, 48] as follows

$$\mathrm{C}{\mathrm{T}}_{\mathrm{\alpha }}=\frac{\mathrm{D}}{{\mathbb{K}}_{\mathrm{\alpha }}}{\int }_{\mathrm{O}}^{\mathrm{\tau }}{\int }_{\mathrm{S}}\left[{\mathrm{T}}_{\mathrm{\alpha }}{\mathrm{q}}^{\ast }-{\mathrm{T}}_{\mathrm{\alpha }}^{\ast }\mathrm{q}\right]\mathrm{dS}\mathrm{d\tau }+\frac{\mathrm{D}}{{\mathbb{K}}_{\mathrm{\alpha }}}{\int }_{\mathrm{O}}^{\mathrm{\tau }}{\int }_{\mathrm{R}}\mathrm{b}{\mathrm{T}}_{\mathrm{\alpha }}^{\ast }\mathrm{dR}\mathrm{d\tau }+{\int }_{\mathrm{R}}{\left.{\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{i}}{\mathrm{T}}_{\mathrm{\alpha }}^{\ast }\right|}_{\mathrm{\tau }=0}$$ E17

which can be expressed in the following form [53].

$$\mathrm{C}{\mathrm{T}}_{\mathrm{\alpha }}={\int }_{\mathrm{S}}\left[{\mathrm{T}}_{\mathrm{\alpha }}{\mathrm{q}}^{\ast }-{\mathrm{T}}_{\mathrm{\alpha }}^{\ast }\mathrm{q}\right]\mathrm{dS}-{\int }_{\mathrm{R}}\frac{{\mathbb{K}}_{\mathrm{\alpha }}}{\mathrm{D}}\frac{\partial {\mathrm{T}}_{\mathrm{\alpha }}^{\ast }}{\mathrm{\partial \tau }}{\mathrm{T}}_{\mathrm{\alpha }}\mathrm{dR}$$ E18

The time derivative of temperature ${\mathrm{T}}_{\mathrm{\alpha }}$ can be approximated as

$$\frac{\partial {\mathrm{T}}_{\mathrm{\alpha }}}{\mathrm{\partial \tau }}\cong \sum _{\mathrm{j}=1}^{\mathrm{N}}{\mathrm{f}}^{\mathrm{j}}{\left(\mathrm{r}\right)}^{\mathrm{j}}{\mathrm{a}}^{\mathrm{j}}\left(\mathrm{\tau }\right).$$ E19

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

Also, we assume that ${\widehat{\mathrm{T}}}_{\mathrm{\alpha }}^{\mathrm{j}}$ is a solution of

$${\nabla }^2{\widehat{\mathrm{T}}}_{\mathrm{\alpha }}^{\mathrm{j}}={\mathrm{f}}^{\mathrm{j}}$$ E20

Thus, Eq. (18) results in the following boundary integral equation [53]

$$\mathrm{C}{\mathrm{T}}_{\mathrm{\alpha }}={\int }_{\mathrm{S}}\left[{\mathrm{T}}_{\mathrm{\alpha }}{\mathrm{q}}^{\ast }-{\mathrm{T}}_{\mathrm{\alpha }}^{\ast }\mathrm{q}\right]\mathrm{dS}+\sum _{\mathrm{j}=1}^{\mathrm{N}}{\mathrm{a}}^{\mathrm{j}}\left(\mathrm{\tau }\right){\mathrm{D}}^{-1}\left(\mathrm{C}{\widehat{\mathrm{T}}}_{\mathrm{\alpha }}^{\mathrm{j}}-{\int }_{\mathrm{S}}\left[{\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{j}}{\mathrm{q}}^{\ast }-{\widehat{\mathrm{q}}}^{\mathrm{j}}{\mathrm{T}}_{\mathrm{\alpha }}^{\ast }\right]\mathrm{dS}\right)$$ E21

where

$${\widehat{\mathrm{q}}}^{\mathrm{j}}=-{\mathbb{K}}_{\mathrm{\alpha }}\frac{\partial {\widehat{\mathrm{T}}}_{\mathrm{\alpha }}^{\mathrm{j}}}{\mathrm{\partial n}}$$ E22

and

$${\mathrm{a}}^{\mathrm{j}}\left(\mathrm{\tau }\right)=\sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{f}}_{\mathrm{ji}}^{-1}\frac{\mathrm{\partial T}\left({\mathrm{r}}_{\mathrm{i}},\mathrm{\tau }\right)}{\mathrm{\partial \tau }}$$ E23

In which the entries of ${\mathrm{f}}_{\mathrm{ji}}^{-1}$ are the coefficients of ${\mathrm{F}}^{-1}$ which described in [34].

$${\left\{\mathrm{F}\right\}}_{\mathrm{ji}}={\mathrm{f}}^{\mathrm{j}}\left({\mathrm{r}}_{\mathrm{i}}\right)$$ E24

The boundary integral discretization scheme has been applied to (21) with the use of (23), we get [53]

$$\mathrm{C}{\dot{\mathrm{T}}}_{\mathrm{\alpha }}+\mathrm{H}{\mathrm{T}}_{\mathrm{\alpha }}=\mathrm{G}\mathrm{Q}$$ E25

where ${\mathrm{T}}_{\mathrm{\alpha }}$ and $\mathrm{Q}$ are temperature, heat flux vectors and internal heat generation vectors, respectively.

The diffusion matrix can be defined as

$$\mathrm{C}=-\left[\mathrm{H}{\widehat{\mathrm{T}}}_{\mathrm{\alpha }}-\mathrm{G}\widehat{\mathrm{Q}}\right]{\mathrm{F}}^{-1}{\mathrm{D}}^{-1}$$ E26

with

$${\left\{\widehat{\mathrm{T}}\right\}}_{\mathrm{ij}}={\widehat{\mathrm{T}}}^{\mathrm{j}}\left({\mathrm{x}}_{\mathrm{i}}\right)$$ E27

$${\left\{\widehat{\mathrm{Q}}\right\}}_{\mathrm{ij}}={\widehat{\mathrm{q}}}^{\mathrm{j}}\left({\mathrm{x}}_{\mathrm{i}}\right)$$ E28

For solving (25) numerically, we interpolate the functions ${\mathrm{T}}_{\mathrm{\alpha }}$ and $\mathrm{q}$ as

$${\mathrm{T}}_{\mathrm{\alpha }}=\left(1-\mathrm{\theta }\right){\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{m}}+\mathrm{\theta }{\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{m}+1}$$ E29

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

where $0\le \mathrm{\theta }=\frac{\mathrm{\tau }-{\mathrm{\tau }}^{\mathrm{m}}}{{\mathrm{\tau }}^{\mathrm{m}+1}-{\mathrm{\tau }}^{\mathrm{m}}}\le 1$ .

The time derivative of (29) can be expressed as

$${\dot{\mathrm{T}}}_{\mathrm{\alpha }}=\frac{\mathrm{d}{\mathrm{T}}_{\mathrm{\alpha }}}{\mathrm{d\theta }}\frac{\mathrm{d\theta }}{\mathrm{d\tau }}=\frac{{\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{m}+1}-{\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{m}}}{{\mathrm{\tau }}^{\mathrm{m}+1}-{\mathrm{\tau }}^{\mathrm{m}}}=\frac{{\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{m}+1}-{\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{m}}}{∆{\mathrm{\tau }}^{\mathrm{m}}}$$ E31

By substituting from Eqs. (29)–(31) into Eq. (25), we obtain

$$\left(\frac{\mathrm{C}}{∆{\mathrm{\tau }}^{\mathrm{m}}}+\mathrm{\theta H}\right){\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{m}+1}-\mathrm{\theta G}{\mathrm{Q}}^{\mathrm{m}+1}=\left(\frac{\mathrm{C}}{∆{\mathrm{\tau }}^{\mathrm{m}}}-\left(1-\mathrm{\theta }\right)\mathrm{H}\right){\mathrm{T}}_{\mathrm{\alpha }}^{\mathrm{m}}+\left(1-\mathrm{\theta }\right)\mathrm{G}{\mathrm{Q}}^{\mathrm{m}}$$ E32

Making use of initial conditions and boundary conditions at $∆{\mathrm{\tau }}^{\mathrm{m}}$ and considering the previous time step solution as initial values for next step, we get

$$\mathbb{a}\mathrm{X}=\mathbb{b}$$ E33

The Adaptive Smoothing and Prolongation Algebraic Multigrid (aSP-AMG) method, which uses an adaptive Factorized Sparse Approximate Inverse (aFSAI) [67] preconditioner as high performance technique that has been implemented efficiently in Matlab (R2018a) for solving the resulting simultaneous linear algebraic systems (33).

3.2 BEM simulation for micropolar thermoviscoelastic fields

According to the weighted residual method, we can write the differential Eqs. (1) and (2) in the following integral form

$${\int }_{\mathrm{R}}\left({\mathrm{\sigma }}_{\mathrm{ij},\mathrm{j}}+{\mathrm{U}}_{\mathrm{i}}\right){\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dR}=0$$ E34

$${\int }_{\mathrm{R}}\left({\mathrm{m}}_{\mathrm{ij},\mathrm{j}}+{\mathrm{\epsilon }}_{\mathrm{ijk}}{\mathrm{\sigma }}_{\mathrm{jk}}+{\mathrm{V}}_{\mathrm{i}}\right){\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dR}=0$$ E35

where

$${\mathrm{U}}_{\mathrm{i}}=\mathrm{\rho }{\mathrm{F}}_{\mathrm{i}}-\mathrm{\rho }{\ddot{\mathrm{u}}}_{\mathrm{i}},{\mathrm{V}}_{\mathrm{i}}=\mathrm{\rho }\left({\mathrm{M}}_{\mathrm{i}}-\mathrm{J}{\ddot{\mathrm{\omega }}}_{\mathrm{i}}\right)$$ E36

The boundary conditions are

$${\mathrm{u}}_{\mathrm{i}}={\overline{\mathrm{u}}}_{\mathrm{i}}\mathrm{on}{\mathrm{S}}_1$$ E37

$${\mathrm{\lambda }}_{\mathrm{i}}={\mathrm{\sigma }}_{\mathrm{ij}}{\mathrm{n}}_{\mathrm{j}}={\overline{\mathrm{\lambda }}}_{\mathrm{i}}\mathrm{on}{\mathrm{S}}_2$$ E38

$${\mathrm{\omega }}_{\mathrm{i}}={\overline{\mathrm{\omega }}}_{\mathrm{i}}\mathrm{on}{\mathrm{S}}_3$$ E39

$${\mathrm{\mu }}_{\mathrm{i}}={\mathrm{m}}_{\mathrm{ij}}{\mathrm{n}}_{\mathrm{j}}={\overline{\mathrm{\mu }}}_{\mathrm{i}}\mathrm{on}{\mathrm{S}}_4$$ E40

By integrating by parts the first term of Eqs. (34) and (35), we obtain

$$-{\int }_{\mathrm{R}}{\mathrm{\sigma }}_{\mathrm{ij}}{\mathrm{u}}_{\mathrm{i},\mathrm{j}}^{\ast }\mathrm{dR}+{\int }_{\mathrm{R}}{\mathrm{U}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dR}=-{\int }_{{\mathrm{S}}_2}{\mathrm{\lambda }}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dS}$$ E41

$$-{\int }_{\mathrm{R}}{\mathrm{m}}_{\mathrm{ij}}{\mathrm{\omega }}_{\mathrm{i},\mathrm{j}}^{\ast }\mathrm{dR}+{\int }_{\mathrm{R}}{\mathrm{\epsilon }}_{\mathrm{ijk}}{\mathrm{\sigma }}_{\mathrm{jk}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dR}+{\int }_{\mathrm{R}}{\mathrm{V}}_{\mathrm{i}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dR}=-{\int }_{{\mathrm{S}}_4}{\mathrm{\mu }}_{\mathrm{i}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dS}$$ E42

On the basis of Huang and Liang [64], we can write

$$\begin{array}{l}-{\int }_{\mathrm{R}}{\mathrm{\sigma }}_{\mathrm{ij},\mathrm{j}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dR}+{\int }_{\mathrm{R}}\left({\mathrm{m}}_{\mathrm{ij},\mathrm{j}}+{\mathrm{\epsilon }}_{\mathrm{ijk}}{\mathrm{\sigma }}_{\mathrm{jk}}\right){\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dR}+{\int }_{\mathrm{R}}{\mathrm{U}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dR}+{\int }_{\mathrm{R}}{\mathrm{V}}_{\mathrm{i}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dR}\\ ={\int }_{{\mathrm{S}}_2}\left({\mathrm{\lambda }}_{\mathrm{i}}-{\overline{\mathrm{\lambda }}}_{\mathrm{i}}\right){\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dS}+{\int }_{{\mathrm{S}}_1}\left({\overline{\mathrm{u}}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{i}}\right){\mathrm{\lambda }}_{\mathrm{i}}^{\ast }\mathrm{dS}+{\int }_{{\mathrm{S}}_4}\left({\mathrm{\mu }}_{\mathrm{i}}-{\overline{\mathrm{\mu }}}_{\mathrm{i}}\right){\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dS}+{\int }_{{\mathrm{S}}_3}\left({\overline{\mathrm{\omega }}}_{\mathrm{i}}-{\mathrm{\omega }}_{\mathrm{i}}\right){\mathrm{\mu }}_{\mathrm{i}}^{\ast }\mathrm{dS}\end{array}$$ E43

By integrating by parts, the left-hand side of (43) can be written as

$$\begin{array}{l}-{\int }_{\mathrm{R}}{\mathrm{\sigma }}_{\mathrm{ij}}{\mathrm{\epsilon }}_{\mathrm{ij}}^{\ast }\mathrm{dR}-{\int }_{\mathrm{R}}{\mathrm{m}}_{\mathrm{ij},\mathrm{j}}{\mathrm{\omega }}_{\mathrm{i},\mathrm{j}}^{\ast }\mathrm{dR}+{\int }_{\mathrm{R}}{\mathrm{U}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dR}+{\int }_{\mathrm{R}}{\mathrm{V}}_{\mathrm{i}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dR}\\ =-{\int }_{{\mathrm{S}}_2}{\overline{\mathrm{\lambda }}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dS}-{\int }_{{\mathrm{S}}_1}{\mathrm{\lambda }}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dS}+{\int }_{{\mathrm{S}}_1}\left({\overline{\mathrm{u}}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{i}}\right){\mathrm{\lambda }}_{\mathrm{i}}^{\ast }\mathrm{dS}-{\int }_{{\mathrm{S}}_4}{\overline{\mathrm{\mu }}}_{\mathrm{i}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dS}-{\int }_{{\mathrm{S}}_3}{\mathrm{\mu }}_{\mathrm{i}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dS}+{\int }_{{\mathrm{S}}_3}\left({\overline{\mathrm{\omega }}}_{\mathrm{i}}-{\mathrm{\omega }}_{\mathrm{i}}\right){\mathrm{\mu }}_{\mathrm{i}}^{\ast }\mathrm{dS}\end{array}$$ E44

According to Eringen [68], the elastic and couple stresses can be written in the following form

$${\mathrm{\sigma }}_{\mathrm{ij}}={\mathbb{A}}_{\text{ijkl}}{\mathrm{\epsilon }}_{\mathrm{kl}},{\mathrm{m}}_{\mathrm{ij}}={\mathbb{B}}_{\text{ijkl}}{\mathrm{\omega }}_{\mathrm{k},\mathrm{l}}$$ E45

where ${\mathbb{A}}_{\text{ijkl}}={\mathbb{A}}_{\text{klij}}\text{and}{\mathbb{B}}_{\text{ijkl}}={\mathbb{B}}_{\text{klij}}$ as shown in [68].

Hence, Eq. (44) can be re-expressed as [53]

$$\begin{array}{l}-{\int }_{\mathrm{R}}{\mathrm{\sigma }}_{\mathrm{ij}}^{\ast }{\mathrm{\epsilon }}_{\mathrm{ij}}\mathrm{dR}-{\int }_{\mathrm{R}}{\mathrm{m}}_{\mathrm{ij},\mathrm{j}}^{\ast }{\mathrm{\omega }}_{\mathrm{i},\mathrm{j}}\mathrm{dR}+{\int }_{\mathrm{R}}{\mathrm{U}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dR}+{\int }_{\mathrm{R}}{\mathrm{V}}_{\mathrm{i}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dR}\\ =-{\int }_{{\mathrm{S}}_2}{\overline{\mathrm{\lambda }}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dS}-{\int }_{{\mathrm{S}}_1}{\mathrm{\lambda }}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{i}}^{\ast }\mathrm{dS}+{\int }_{{\mathrm{S}}_1}\left({\overline{\mathrm{u}}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{i}}\right){\mathrm{\lambda }}_{\mathrm{i}}^{\ast }\mathrm{dS}-{\int }_{{\mathrm{S}}_4}{\overline{\mathrm{\mu }}}_{\mathrm{i}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dS}-{\int }_{{\mathrm{S}}_3}{\mathrm{\mu }}_{\mathrm{i}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }\mathrm{dS}+{\int }_{{\mathrm{S}}_3}\left({\overline{\mathrm{\omega }}}_{\mathrm{i}}-{\mathrm{\omega }}_{\mathrm{i}}\right){\mathrm{\mu }}_{\mathrm{i}}^{\ast }\mathrm{dS}\end{array}$$ E46

By applying integration by parts again, the left-hand side of (46) can be written as [53]

$${\int }_{\mathrm{R}}{\mathrm{\sigma }}_{\mathrm{ij},\mathrm{j}}^{\ast }{\mathrm{u}}_{\mathrm{i}}\mathrm{dR}+{\int }_{\mathrm{R}}\left({\mathrm{m}}_{\mathrm{ij},\mathrm{j}}^{\ast }+{\mathrm{\epsilon }}_{\mathrm{ijk}}{\mathrm{\sigma }}_{\mathrm{jk}}^{\ast }\right){\mathrm{\omega }}_{\mathrm{i}}\mathrm{dR}=-{\int }_{\mathrm{S}}{\mathrm{u}}_{\mathrm{i}}^{\ast }{\mathrm{\lambda }}_{\mathrm{i}}\mathrm{dS}-{\int }_{\mathrm{S}}{\mathrm{\omega }}_{\mathrm{i}}^{\ast }{\mathrm{\mu }}_{\mathrm{i}}\mathrm{dS}+{\int }_{\mathrm{S}}{\mathrm{\lambda }}_{\mathrm{i}}^{\ast }{\mathrm{u}}_{\mathrm{i}}\mathrm{dS}+{\int }_{\mathrm{S}}{\mathrm{\mu }}_{\mathrm{i}}^{\ast }{\mathrm{\omega }}_{\mathrm{i}}\mathrm{dS}$$ E47

The obtained weighting functions for ${\mathrm{U}}_{\mathrm{i}}={∆}^{\mathrm{n}}$ and ${\mathrm{V}}_{\mathrm{i}}=0$ along ${\mathrm{e}}_{\mathrm{l}}$ were first used as follows:

$${\mathrm{\sigma }}_{\mathrm{lj},\mathrm{j}}^{\ast }+{∆}^{\mathrm{n}}{\mathrm{e}}_{\mathrm{l}}=0$$ E48

$${\mathrm{m}}_{\mathrm{ij},\mathrm{j}}^{\ast }+{\mathrm{\epsilon }}_{\mathrm{ijk}}{\mathrm{\sigma }}_{\mathrm{jk}}^{\ast }=0$$ E49

According to Dragos [69], the fundamental solutions can be written as

$${\mathrm{u}}_{\mathrm{i}}^{\ast }={\mathrm{u}}_{\mathrm{li}}^{\ast }{\mathrm{e}}_{\mathrm{l}},{\mathrm{\omega }}_{\mathrm{i}}^{\ast }={\mathrm{\omega }}_{\mathrm{li}}^{\ast }{\mathrm{e}}_{\mathrm{l}},{\mathrm{\lambda }}_{\mathrm{i}}^{\ast }={\mathrm{\lambda }}_{\mathrm{li}}^{\ast }{\mathrm{e}}_{\mathrm{l}},{\mathrm{\mu }}_{\mathrm{i}}^{\ast }={\mathrm{\mu }}_{\mathrm{li}}^{\ast }{\mathrm{e}}_{\mathrm{l}},$$ E50

The obtained weighting functions for ${\mathrm{U}}_{\mathrm{i}}=0$ and ${\mathrm{V}}_{\mathrm{i}}={∆}^{\mathrm{n}}$ along ${\mathrm{e}}_{\mathrm{l}}$ were next used as follows:

$${\mathrm{\sigma }}_{\mathrm{ij},\mathrm{j}}^{\ast \ast }=0$$ E51

$${\mathrm{m}}_{\mathrm{lj},\mathrm{j}}^{\ast \ast }+{\mathrm{\epsilon }}_{\mathrm{ljk}}{\mathrm{\sigma }}_{\mathrm{jk}}^{\ast \ast }+{∆}^{\mathrm{n}}{\mathrm{e}}_{\mathrm{l}}=0$$ E52

The fundamental solutions that have been obtained analytically by Dragos [69] can be written as

$${\mathrm{u}}_{\mathrm{i}}^{\ast }={\mathrm{u}}_{\mathrm{li}}^{\ast \ast }{\mathrm{e}}_{\mathrm{l}},{\mathrm{\omega }}_{\mathrm{i}}^{\ast }={\mathrm{\omega }}_{\mathrm{li}}^{\ast \ast }{\mathrm{e}}_{\mathrm{l}},{\mathrm{\lambda }}_{\mathrm{i}}^{\ast }={\mathrm{\lambda }}_{\mathrm{li}}^{\ast \ast }{\mathrm{e}}_{\mathrm{l}},{\mathrm{\mu }}_{\mathrm{i}}^{\ast }={\mathrm{\mu }}_{\mathrm{li}}^{\ast \ast }{\mathrm{e}}_{\mathrm{l}},$$ E53

Using the above two sets of weighting functions into (47) we have

$${\mathrm{C}}_{\mathrm{li}}^{\mathrm{n}}{\mathrm{u}}_{\mathrm{i}}^{\mathrm{n}}=-{\int }_{\mathrm{S}}{\mathrm{\lambda }}_{\mathrm{li}}^{\ast }{\mathrm{u}}_{\mathrm{i}}\mathrm{dS}-{\int }_{\mathrm{S}}{\mathrm{\mu }}_{\mathrm{li}}^{\ast }{\mathrm{\omega }}_{\mathrm{i}}\mathrm{dS}+{\int }_{\mathrm{S}}{\mathrm{u}}_{\mathrm{li}}^{\ast }{\mathrm{\lambda }}_{\mathrm{i}}\mathrm{dS}+{\int }_{\mathrm{S}}{\mathrm{\omega }}_{\mathrm{li}}^{\ast }{\mathrm{\mu }}_{\mathrm{i}}\mathrm{dS}$$ E54

$${\mathrm{C}}_{\mathrm{li}}^{\mathrm{n}}{\mathrm{\omega }}_{\mathrm{i}}^{\mathrm{n}}=-{\int }_{\mathrm{S}}{\mathrm{\lambda }}_{\mathrm{li}}^{\ast \ast }{\mathrm{u}}_{\mathrm{i}}\mathrm{dS}-{\int }_{\mathrm{S}}{\mathrm{\mu }}_{\mathrm{li}}^{\ast \ast }{\mathrm{\omega }}_{\mathrm{i}}\mathrm{dS}+{\int }_{\mathrm{S}}{\mathrm{u}}_{\mathrm{li}}^{\ast \ast }{\mathrm{\lambda }}_{\mathrm{i}}\mathrm{dS}+{\int }_{\mathrm{S}}{\mathrm{\omega }}_{\mathrm{li}}^{\ast \ast }{\mathrm{\mu }}_{\mathrm{i}}\mathrm{dS}$$ E55

Thus, we can write

$${\mathrm{C}}^{\mathrm{n}}{\mathbb{q}}^{\mathrm{n}}=-{\int }_{\mathrm{S}}{\mathbb{p}}^{\ast }\mathbb{q}\mathrm{dS}+{\int }_{\mathrm{S}}{\mathbb{q}}^{\ast }\mathbb{p}\mathrm{dS}$$ E56

where

$${\mathit{C}}^{\mathit{n}}=\left[\begin{array}{ccc}{\mathit{C}}_{11}& {\mathit{C}}_{12}& 0\\ {\mathit{C}}_{21}& {\mathit{C}}_{22}& 0\\ 0& 0& 0\end{array}\right], {\mathbb{q}}^{\ast }=\left[\begin{array}{ccc}{\mathit{u}}_{11}^{\ast }& {\mathit{u}}_{12}^{\ast }& {\mathit{\omega }}_{13}^{\ast }\\ {\mathit{u}}_{21}^{\ast }& {\mathit{u}}_{22}^{\ast }& {\mathit{\omega }}_{23}^{\ast }\\ {\mathit{u}}_{31}^{\ast \ast }& {\mathit{u}}_{32}^{\ast \ast }& {\mathit{\omega }}_{33}^{\ast \ast }\end{array}\right], {\mathbb{p}}^{\ast }=\left[\begin{array}{ccc}{\mathit{\lambda }}_{11}^{\ast }& {\mathit{\lambda }}_{12}^{\ast }& {\mathit{\mu }}_{13}^{\ast }\\ {\mathit{\lambda }}_{21}^{\ast }& {\mathit{\lambda }}_{22}^{\ast }& {\mathit{\mu }}_{23}^{\ast }\\ {\mathit{\lambda }}_{31}^{\ast \ast }& {\mathit{\lambda }}_{32}^{\ast \ast }& {\mathit{\mu }}_{33}^{\ast \ast }\end{array}\right],\mathbb{q}=\left[\begin{array}{c}{\mathit{u}}_1\\ {\mathit{u}}_2\\ {\mathit{\omega }}_3\end{array}\right],\mathbb{p}=\left[\begin{array}{c}{\mathit{\lambda }}_1\\ {\mathit{\lambda }}_2\\ {\mathit{\mu }}_3\end{array}\right]$$

In order to solve (56) numerically, we construct the following functions

$$\mathrm{q}=\mathrm{\psi }{\mathrm{q}}^{\mathrm{j}},\mathrm{p}=\mathrm{\psi }{\mathrm{p}}^{\mathrm{j}}$$ E57

substituting above functions into (56) and discretizing the boundary, we obtain

$${\mathrm{C}}^{\mathrm{n}}{\mathbb{q}}^{\mathrm{n}}=\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{e}}}\left[-{\int }_{{\mathrm{\Gamma }}_{\mathrm{j}}}{\mathbb{p}}^{\ast }\mathrm{\psi }\mathrm{d\Gamma }\right]{\mathbb{q}}^{\mathrm{j}}+\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{e}}}\left[{\int }_{{\mathrm{\Gamma }}_{\mathrm{j}}}{\mathbb{q}}^{\ast }\mathrm{\psi }\mathrm{d\Gamma }\right]{\mathbb{p}}^{\mathrm{j}}$$ E58

Equation after integration may be expressed as

$${\mathrm{C}}^{\mathrm{i}}{\mathbb{q}}^{\mathrm{i}}=-\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{e}}}{\widehat{\mathbb{H}}}^{\mathrm{ij}}{\mathbb{q}}^{\mathrm{j}}+\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{e}}}{\widehat{\mathbb{G}}}^{\mathrm{ij}}{\mathbb{p}}^{\mathrm{j}}$$ E59

which can be expressed as

$$\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{e}}}{\mathbb{H}}^{\mathrm{ij}}{\mathbb{q}}^{\mathrm{j}}=\sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{e}}}{\widehat{\mathbb{G}}}^{\mathrm{ij}}{\mathbb{p}}^{\mathrm{j}}$$ E60

where

$${\mathbb{H}}^{\mathrm{ij}}=\left\{\begin{array}{l}{\widehat{\mathbb{H}}}^{\mathrm{ij}}\text{if}\mathrm{i}\ne \mathrm{j}\\ {\widehat{\mathbb{H}}}^{\mathrm{ij}}+{\mathrm{C}}^{\mathrm{i}}\text{if}\mathrm{i}=\mathrm{j}\end{array}\right.$$ E61

Thus, we can write the following system of matrix equation as

$$\mathbb{HQ}=\mathbb{GP}$$ E62

Hence, we get the following system of linear algebraic equations

$$\mathbb{A}\mathbb{X}=\mathbb{B}$$ E63

4. A new boundary element technique for simulation and optimization of solid deformable bodies under different loads

In order to solve (63), we apply adaptive smoothing and prolongation algebraic multigrid (aSP-AMG) based on adaptive Factorized Sparse Approximate Inverse (aFSAI) as described in [67] for solving the resulting simultaneous linear algebraic system (63) in Matlab (R2018a).

B-spline basis functions are used as weights in the same manner as Bézier basis functions. Spline curves can be expressed in terms of $\mathrm{k}-\text{order}$ B-spline basis function. All B-spline basis functions are assumed to have their domain on [0,1]. B-splines basis functions are a more general type of curve than Bezier curves, where each control point ${\mathrm{P}}_{\mathrm{i}}$ of $\mathrm{i}+1$ control points $\left({\mathrm{P}}_0,{\mathrm{P}}_1,{\mathrm{P}}_2,\dots ,{\mathrm{P}}_{\mathrm{i}}\right)$ is connected with a basis function ${\mathrm{N}}_{\mathrm{i},\mathrm{k}}$ , the knots are the points that subdivide the domain [0,1] into knot spans. Also, each B-spline basis function is non-zero on the entire interval.

The efficiency of our numerical modeling technique has been improved using a nonuniform rational B-spline curve (NURBS) to decrease the computation time and model optimized boundary where it reduces the number of control points and provides the flexibility to design a large variety of shapes.

The considered NURBS can be defined as follows

$$\mathrm{C}\left(\mathrm{t}\right)=\frac{{\sum }_{\mathrm{i}=0}^{\mathrm{n}}{\mathrm{N}}_{\mathrm{i},\mathrm{o}}\left(\mathrm{t}\right){\mathrm{\varpi }}_{\mathrm{i}}{\mathrm{P}}_{\mathrm{i}}}{{\sum }_{\mathrm{i}=0}^{\mathrm{n}}{\mathrm{N}}_{\mathrm{i},\mathrm{o}}\left(\mathrm{t}\right){\mathrm{\varpi }}_{\mathrm{i}}}$$ E64

where ${\mathrm{N}}_{\mathrm{i},\mathrm{o}}\left(\mathrm{t}\right)$ and ${\mathrm{\varpi }}_{\mathrm{i}}$ are the B-spline basis functions of order o and the weights of control points ${\mathrm{P}}_{\mathrm{i}}$ , respectively.

The genetic algorithm greatly reduces computing time and computer memory of achieving an optimum solution, so, it can be used for solving multi-objective problems without needing to calculate the sensitivities. The profiles of the considered objects are represented based on the free form deformation (FFD) technique, where the FFD control points are considered as the genes and then the profiles of chromosomes are defined by the sequence of genes. The population is constructed by many individuals (chromosomes), where the fitness functions are evaluated by using the BEM.

Two criteria can be implemented during shape optimization of the solid bodies [70]

I. The minimum global compliance:

$$\mathcal{F}=\frac{1}{2}\underset{\mathrm{S}}{\int }\left(\mathrm{\lambda }·\mathrm{u}\right)\mathrm{dS}$$ E65

II. The minimum boundary equivalent stresses

$$\mathcal{F}=\underset{\mathrm{S}}{\int }{\left(\frac{{\mathrm{\sigma }}_{\mathrm{ij}}}{{\mathrm{\sigma }}_0}\right)}^{\mathrm{n}}\mathrm{dS}$$ E66

${\mathrm{\sigma }}_{\mathrm{ij}},{\mathrm{\sigma }}_0\text{and}\mathrm{n}$ are equivalent boundary stresses, reference stress and natural number, respectively, where the greater value of $\mathrm{n}$ increases the speed of convergence of the functional (66). By minimizing the functional (66) ${\mathrm{\sigma }}_{i\mathrm{j}}$ are closer to ${\mathrm{\sigma }}_0$ .

In order to find the optimal boundary conditions for temperature the following functional can be applied

$$\mathcal{F}=\underset{\mathrm{S}}{\int }{\left(\frac{\mathrm{u}}{{\mathrm{u}}_0}\right)}^{\mathrm{n}}\mathrm{dS}$$ E67

where $u$ and $u_0$ are boundary displacement and reference displacement, respectively. Minimization of the functional (67) reduces displacements on the selected part of the boundary.

In order to identify unknown inner boundary, we use the following functional

$$\mathcal{F}=\mathrm{\delta }\sum _{\mathrm{k}=1}^{\mathrm{M}}\left({\mathrm{u}}^{\mathrm{k}}-{\widehat{\mathrm{u}}}^{\mathrm{k}}\right)+\mathrm{\eta }\sum _{\mathrm{l}=1}^{\mathrm{N}}\left({\mathrm{T}}^{\mathrm{l}}-{\widehat{\mathrm{T}}}^{\mathrm{l}}\right)$$ E68

where ${\widehat{\mathrm{u}}}^{\mathrm{k}}$ and ${\widehat{\mathrm{T}}}^{\mathrm{l}}$ are measured displacements and temperatures in boundary points $k$ and $l$ respectively, ${\mathrm{u}}^{\mathrm{k}}$ and ${\mathrm{T}}^{\mathrm{l}}$ are computed displacements and temperatures in boundary points $\mathrm{k}$ and $\mathrm{l}$ respectively, $\mathrm{\delta }$ and $\mathrm{\eta }$ are weight coefficients, and $\mathrm{M}$ and $\mathrm{N}$ are numbers of sensors.

5. Numerical examples and discussion of results

For illustration of the theoretical results of our proposed technique from the preceding sections, two numerical examples are analyzed below. The first example is the cantilever beam with inferior corner load, the second is the Michell-type structure, where the material has the following physical data [58]:

The elasticity tensor

$$C_{\mathit{\text{ijkl}}}=\left(\begin{array}{cccccc}60.23& 18.67& 18.96& -7.69& 15.60& -25.28\\ 18.67& 21.26& 9.36& -3.74& 4.21& -8.47\\ 18.96& 9.36& 47.04& -8.82& 15.28& -8.31\\ -7.69& -3.74& -8.82& 10.18& -9.54& 5.69\\ 15.60& 4.21& 15.28& -9.54& 21.19& -8.54\\ -25.28& -8.47& -8.31& 5.69& -8.54& 20.75\end{array}\right)$$ E69

$\mathrm{p}=25\mathrm{MPa}$ , and $\mathrm{∆t}=0.0006\mathrm{s}$ .

Example 1. Cantilever beam structure.

As a practical example, the proposed algorithm is applied on the cantilever beam with inferior corner load $\mathrm{P}=100\mathrm{N}/\mathrm{mm}$ . The geometry of the cantilever beam is illustrated in Figure 1. The initial boundary mesh of the cantilever beam composed of 14 quadratic elements is also illustrated in Figure 2. The BEM grid is composed of 76 nodes along $x$ direction and 51 nodes along $y$ direction. These mesh parameters were obtained after convergence analysis. In the process of optimization, the cantilever beam structure optimization results are presented in Figure 3 from initial to final structure for different iterations.

Figure 1.

Cantilever beam structure geometry.

Figure 2.

Initial boundary of the cantilever beam structure.

Figure 3.

Cantilever beam optimization process from initial to final structure for different iterations.

The present measured boundary element method (BEM) optimization results of the first example are compared in Figure 4 with measured finite difference method (FDM) optimization results obtained by Itzá et al. [71] and measured finite element method (FEM) optimization results obtained using the software package COMSOL Multiphysics, version 5.4. It is clear from this figure that the BEM results obtained by the proposed technique are in excellent agreement with the FDM results [71] and FEM results of the COMSOL Multiphysics.

Figure 4.

Final cantilever beam structure for BEM, FDM and FEM.

Table 1 shows that our proposed BEM modeling of cantilever beam with inferior corner load drastically reduces the manpower needed for modeling and computer resources needed for the calculation in comparison with the calculated results based on the FDM and FEM.

Table 1.

Comparison of computer resources needed for FDM, FEM and BEM modeling of cantilever beam structure.

Example 2. Michell-type structure.

As application example, we use a beam with a mid-span load $\left(\mathrm{P}=100\mathrm{N}/\mathrm{mm}\right)$ (Michell-type structure) as shown in Figure 5. The initial boundary mesh of the Michell-type structure composed of 40 quadratic discontinuous elements is also illustrated in Figure 6. The BEM grid is composed of 76 nodes along $x$ direction, and 51 nodes along $y$ direction. This grid density was obtained after convergence analysis.

Figure 5.

Michel-type structure geometry.

Figure 6.

Initial boundary of the Michel-type structure.

Figure 7 shows the cantilever beam optimization process from initial to final structure for different iterations.

Figure 7.

Michell-type structure optimization process from initial to final structure for different iterations.

The optimization results of the second example obtained with the proposed BEM are compared in Figure 8 with FDM optimization results [71] and FEM optimization results of COMSOL Multiphysics software, version 5.4. It is clear from this figure that our BEM results obtained by the proposed technique are in excellent agreement with the FDM and FEM results.

Figure 8.

Final Michell-type structure for BEM, FDM and FEM.

Table 2 shows that our proposed BEM modeling of Michell-type structure dramatically reduces the computer resources necessary to calculate our proposed modeling in comparison with the calculated results based on the FDM and FEM.

Table 2.

Comparison of computer resources needed for FDM, FEM and BEM modeling of Michell-type structure.

6. Conclusion

In the present paper, we propose a new theory called nonlinear micropolar thermoviscoelasticity involving three temperatures. A new mathematical modeling of nonlinear generalized micropolar thermoviscoelasticity problem. A new boundary element technique for simulation and optimization problems of mechanics of solid deformable bodies is implemented based on genetic algorithm (GA), free form deformation (FFD) method and nonuniform rational B-spline curve (NURBS) as the global optimization technique for solving complex simulation and optimization problems associated with the proposed theory. FFD is an efficient and versatile parameterization technique for treating shape optimization problems with complex shapes. It is implemented for simulation and optimization of the shape. In the formulation of the considered problem, solutions are obtained for specific arbitrary parameters which are the control points positions in the considered problem, the profiles of the considered objects are determined by FFD method, where the FFD control points positions are treated as genes, and then the chromosomes profiles are defined with the genes sequence. The population is founded by a number of individuals (chromosomes), where the objective functions of individuals are determined by the boundary element method (BEM). Due to the large amount of computer resources required by the FDM and FEM, our proposed BEM model can be applied to a wide range of simulation and optimization problems related with our proposed theory. The numerical results demonstrate the validity, accuracy and efficiency of our proposed technique.

α , α ¯ , α ̿ , α ˇ	micro-polar constants
β ij	stress–temperature coefficients
δ ij	Kronecker delta i j = 1 2
δ , η	weight coefficients
ε ij	strain tensor
ε ijk	alternate tensor
ϵ ij	micro-strain tensor
λ	tractions
ϑ 0	viscoelastic relaxation time
ϖ	weights of control points
ρ	material density
σ ij	force stress tensor
σ 0	reference stress
τ	time
τ 0 , τ 1 , τ 2	relaxation times
ω i	micro-rotation vector
ℵ	= 1 + ϑ 0 ∂ ∂ τ viscoelastic constant
b	internal heat generation vector
c	specific heat capacity
C ijkl	constant elastic moduli
e	= ε kk = ϵ kk dilatation
e lij	piezoelectric tensor
F i	mass force vector
J	micro-inertia coefficient
J	current density vector
K e	electron conductive coefficients
K i	ion conductive coefficients
K r	phonon conductive coefficients
k , l	boundary points
M , N	sensors numbers
M i	mass couple vector
m ij	couple stress tensor
P	total energy of unit mass
P e	= c e T e Electron energy
P i	= c i T i Ion energy
P r	= 1 ρ c r T r 4 Phonon energy
ℙ	values vector of tractions and couple stress
p i	singular points
Q	heat flux vectors
ℚ	values vector of displacements and microrotations
R	problem’s boundary
S	problem’s domain
T e	electron temperature
T i	ion temperature
T r	phonon temperature
u	boundary displacement
u i	displacement vector
u 0	reference displacement
u k	computed displacements
u ̂ k	measured displacements
W ei	electron-ion energy coefficient
W er	electron–phonon energy coefficient