A new boundary element algorithm for modeling and simulation of nonlinear thermal stresses in micropolar FGA composites with temperature-dependent properties

The main aim of this article is to develop a new boundary element method (BEM) algorithm to model and simulate the nonlinear thermal stresses problems in micropolar functionally graded anisotropic (FGA) composites with temperature-dependent properties. Some inside points are chosen to treat the nonlinear terms and domain integrals. An integral formulation which is based on the use of Kirchhoff transformation is firstly used to simplify the transient heat conduction governing equation. Then, the residual nonlinear terms are carried out within the current formulation. The domain integrals can be effectively treated by applying the Cartesian transformation method (CTM). In the proposed BEM technique, the nonlinear temperature is computed on the boundary and some inside domain integral. Then, nonlinear displacement can be calculated at each time step. With the calculated temperature and displacement distributions, we can obtain the values of nonlinear thermal stresses. The efficiency of our proposed methodology has been improved by using the communication-avoiding versions of the Arnoldi (CA-Arnoldi) preconditioner for solving the resulting linear systems arising from the BEM to reduce the iterations number and computation time. The numerical outcomes establish the influence of temperature-dependent properties on the nonlinear temperature distribution, and investigate the effect of the functionally graded parameter on the nonlinear displacements and thermal stresses, through the micropolar FGA composites with temperature-dependent properties. These numerical outcomes also confirm the validity, precision and effectiveness of the proposed modeling and simulation methodology.

tailored in accordance with the mechanical-technological properties required in different working conditions. FGMs are excellent advanced materials that change the manufacturing world for the better. FGMs have received a large amount of attention due to its ability to produce materials with tailored properties which are suitable candidates for several hightech applications such as graded structures on the atomic scale, graded hip implants, structural walls, sports equipment, design attractive interference colors for automobiles, etc. Then, the number of publications in this area of research has been growing exponentially in the past 20 years [5][6][7][8][9][10].
Temperature-dependent materials are those substances where the material properties depending upon the temperature at which they are measured, such as ice, bitumen, asphalt concrete, titanium, stainless steel, silver, gold, platinum, granular materials and reverse solubility polymers which are soluble in cold water and insoluble in hot water [11][12][13][14][15][16][17].
In the present paper, a novel boundary element algorithm is proposed for modeling and simulation of nonlinear thermal stresses problems in micropolar FGA composites with temperature-dependent properties. In the proposed BEM formulation, the residual nonlinear terms are treated by using the Kirchhoff transformation and the domain integrals are treated by using the CTM. Then, at each time step the nonlinear temperature, displacements and thermal stresses are calculated at boundary nodes, and a few internal points which are considered to be utilized as initial values for another time step. The numerical findings investigate the effects of temperature-dependent properties and functionally graded parameter. These numerical findings also confirm the validity, precision and effectiveness of the proposed modeling and simulation methodology.

Formulation of the problem
Geometry of the current computational problem is shown in Fig. 1. The governing equations for nonlinear thermal stresses in micropolar composites with temperature-dependent properties can be expressed as ∂�(X, t) ∂t Now, the right-hand side of (8) can be decomposed into linear and nonlinear parts as where 0 , ρ 0 and c 0 are thermal conductivity, density and specific heat, respectively, at T 0 . where the nonlinear term can be written as According to [15,16], Eq. (9) can be written as where The integral equation which corresponds to Eq. (11) can be expressed as The fundamental solution and its normal derivative, respectively, can be written as where a 0 = 0 ρ 0 c 0 and H is the Heaviside function. The time integrals which corresponds to (14) and (15) can be computed analytically as t n+1 t n � * (P, t n+1 ; Q, τ)dτ = 1 4π a 0 Ei r 2 4a 0 �t where the exponential integral function Ei() can be defined as The first domain integral of Eq. (13) contains the nonlinear term, that is or By substituting the midpoint value of h Nl and finite difference expression of ˙ , we can write where Now, we implement the CTM without domain discretisation to evaluate the domain integrals of (23). Thus, the unknown values at M ′ boundary nodes can be computed directly from the following system of matrix equations (17) t n+1 t n q * (P, t n+1 ; Q, τ)dτ = −1 2π a 0 r exp −r 2 4a 0 �t ∂r ∂n where � Ŵ and Q Ŵ are M ′ dimension vectors contain boundary nodal values and q , F a vector depends on previous time step, F NI is a nonlinear term vector depends on unknown internal values, H and G are M ′ × M ′ dimension coefficient matrices. Also, the unknown values at M ′′ internal points may be calculated from the following system of matrix equations where H , G , Ĥ and Ĝ can be computed for all time steps. Also, F , F NI ,F andF NI can be computed using CTM for all time steps.
The CTM method can be implemented to transform several domain integrals into boundary ones [47]. Now, we consider the following two-dimensional regular domain integral By implementing Green's theorem as Now, we can write where Since the integral in (29) cannot be determined analytically, so, we evaluate it numerically by the following integral equation According to Khosravifard and Hematiyan [49], and using (30), the domain integral of (26) can be expressed as where x 1min and x 1max are minimum x 1 and maximum x 1 values, respectively. The composite Gaussian quadrature method is applied to (26) yields Equation (33) can be expressed as where J k and J l are the transformation Jacobian for the kth interval lth interval, respectively, K is the boundary elements number, N and J are the Gaussian integration points numbers of (33) for the outer integral and inner integral, respectively, w i and w j are Gauss points weights.
If p described over a domain-boundary grid with irregularly spaced data. Then, by using the radial point interpolation method (RPIM) [50], the approximation using two-dimensional interpolation of the function p may be written as in which (x 1 , x 2 ) is any arbitrary point, p i is the p value at i and φ i its shape function, M(total number) = M ′ boundary nodes number +M ′′ internal grid points number . In the considered RPIM, the consistent shape functions are constructed using the radial basis functions. According to [50], the function p(x 1 , x 2 ) can be approximated as The considered method is very simple for computation of regular and weakly singular domain integrals because all computations are performed in universal Cartesian coordinates, where kernels are defined by irregularly spaced data.
In order to create the RPIM shape functions, we apply the following Gaussian radial basis functions (GRBFs) where ψ i are radial basis functions (RBFs), n is the RBFs number, m is the polynomial basis functions number and u j (x) , the augmented monomials and α i and b j are unknown coefficients which can be evaluated from the following n linear system of equations.
. . , n and the following m linear constraints From Eqs. (38) and (39), we can write α i and b j in the following form Based on [50], and using (40), Eq. (36) may be expressed as follows where the matrix P is location-and geometry-dependent of boundary nodes and internal nodes.
where φ is the RPIM shape functions vector Now, Eq. (42) can be expressed as in which γ is the geometry-and location-dependent weight vector of grid points and p includes the values p at boundary nodes and internal points.

Regularization of BEM formulations and evaluation of the domain integrals
In the BEM formulation of transient nonlinear thermal stresses problems, there are several regular and singular domain integrals with different kernels should be calculated with boundary-only discretization. Now, we consider the following domain integrals which occur in the integral Eq. (23) where the weakly singular exponential integral function Ei() in (44) can be written as in which   where r i = (x i ) Q − (x i ) P and D 1 (P) has the same value at each iteration of each time step.
The domain integral in (45)

BEM solution for displacement field
The partial differential Eqs. (1) and (2), may be transformed into the following integral equations where (55) in which u * i and ω * i are weighting functions. In the current paper, we considered the following boundary conditions Applying integration by parts to Eqs. (61) and (62), we get Based on Fahmy [41], we can write By using integration by parts for the left-hand side of (71), we obtain On the basis of Fahmy [41], elastic and couple stresses can be written as Hence, Eq. (72) can be rewritten as By applying integration by parts to the left-hand side of (74) for the second time, we get The weighting functions of U i = n and V i = 0 along e l are obtained as follows: As Fahmy put it [41], the fundamental solution may be expressed as Also, the weighting functions of U i = 0 and V i = n along e l can be represented as: Based on Fahmy [41], the fundamental solution can be expressed in the following form Now, by considering the above two sets of weight functions into (75) we obtain Thus, we can write In order to solve (84) numerically, we use the following definitions Substituting the above definitions and properties into (84) and discretizing the boundary into N e elements yield which can be expressed, after integration, as follows: Taking into consideration the following definition Therefore, depends upon the previous definition, Eq. (87) has the following form which may be written as Thus, we obtain A time-stepping algorithm of Fahmy [54,55] based on communication-avoiding Arnoldi preconditioner is implemented in order to obtain the temperature and displacement fields.

Numerical results and discussion
In the current paper, we considered the temperature-dependent pure copper material with the following physical data [61]: The thermal conductivity is The temperature dependent specific heat and density are shown in Tables 1 and 2, respectively.
The proposed BEM technique implemented in the current paper should be applicable to a wide range of nonlinear thermal stresses problems in micropolar FGA structures with temperature-dependent properties.
In the BEM modelling of the considered problem, the boundary has been discretized using 84 linear boundary elements and 404 internal points as shown in Fig. 2.
The effect of functionally graded parameter on the time distributions of the nonlinear displacements and thermal stresses plays an important role during the modeling process. Figure 3 shows the time distribution of the nonlinear temperature for two different cases called temperature dependent (TD) and temperature independent (TID) = 400 1 − T 6000  Table 2 Temperature dependent density of pure cupper [61] T   functionally graded parameter has a tremendous impact on the displacements through the micropolar structures with temperature-dependent properties. Figures 6, 7 and 8 show the time distributions of the nonlinear thermal stresses σ 11 , σ 12 and σ 22 for homogeneous (m = 0) and functionally graded (m = 0.3, 0.6 and 0.9) micropolar structures with temperature-dependent properties. It can be shown from these figures that the functionally graded parameter has a considerable effect on the nonlinear thermal stresses through the micropolar structures with temperature-dependent properties.
The main characteristics of BEM over FDM or FEM, which make it the most appropriate numerical method, for dealing with the considered problem that has been performed are as follows [38,41]: • BEM does not need the internal domain to be discretized. But both FDM and FEM require the discretization of the whole domain. Therefore, BEM is more efficient and easy to use than FDM or FEM. • BEM integration is a smoothing operation and is more numerically stable than differentiation operation of FDM or FEM. Therefore, all variables of the considered problem at any point are more precise than FDM or FEM. • For the solution of the considered open boundary problem, BEM programmers deal with real geometry boundaries, while, FDM or FEM programmers deal with artificial boundaries, which are far away from the real problem. In addition, it is difficult for them to deal with the considered problem.
The efficiency of our proposed technique has been developed by using the communication-avoiding versions of the Arnoldi (CA-Arnoldi) preconditioner for solving the resulting linear systems arising from the BEM to reduce the iterations number and computation time.
The communication-avoiding versions of the Arnoldi (CA-Arnoldi) of Hoemmen [62] that also implemented by Fahmy [41], adaptive smoothing and prolongation algebraic multigrid (aSP-AMG) of Magri et al. [63] that also applied by Fahmy [43] and the regularized of Badahmane [64] which is also used by Fahmy [65] were compared with each other in Table 3. This table reports the iteration number (IT), CPU time, relative residual (RES) and error (ERR) of the tested iteration methods with respect to different values of In the considered special case, the boundary element model of the considered example, the boundary has been discretized using 84 linear boundary elements and 404 internal points as shown in Fig. 2, and the results of temperature and displacements are plotted in Figs. 9, 10, and 11. It can be seen from these figures that the BEM outcomes are in very good agreement with the FDM outcomes of Awrejcewicz and Krysko [66] and FEM of Shakeriaski and Ghodrat [67].

Example 2 Nonlinear thermal stresses in a thick hollow cylinder
In the boundary element model of the considered example, the boundary has been discretized using 44 boundary elements and without internal points as shown in Fig. 12, and the results of nonlinear thermal stresses are plotted in Fig. 13. It can be seen from this figure that the BEM outcomes are in very good agreement with the FDM outcomes of Ahmad et al. [68] and FEM outcomes of Ibrahim and Gadisa [69].

Conclusion
A new boundary element modeling and simulation algorithm is developed based on time-dependent fundamental solutions in order to calculate the nonlinear thermal stresses in micropolar FGA composites, where material properties such as thermal conductivity, density and specific heat are assumed to be temperature-dependent. In the proposed BEM formulation, the boundary is subdivided into boundary elements, and the domain must be subdivided into internal cells, without any connectivity to increase the accuracy of the computation. There is no need to specify particular solutions for evaluating domain integrals, where the Galerkin meshfree method can be used to approximate the integral kernels of the domain integrals. The domain integrals are assessed effectively by the CTM integration weights which are chosen to be constant over each time step. The temperatures at the boundary nodes are the unknowns of the resulting equations system, whereas the internal points temperatures are calculated over all time steps without solving the whole equations system. The proposed iterative algorithm converges quickly and provides highly accurate numerical solution of the nonlinear without need of complex calculation. Comparing the results of the BEM with the FDM and FEM shows that BEM is effective and convenient. Also, it is predicated that BEM can be used in the analysis of strongly nonlinear thermoelastic problems accurately. During our treatment of the problem under consideration, we applied CA-Arnoldi, aSP-AMG and regularized preconditioners. The obtained numerical results demonstrate that CA-Arnoldi preconditioner has better performance than aSP-AMG and regularized preconditioners. The computational validity, accuracy and efficiency of the proposed technique were demonstrated. Hence, this paper demonstrates that, the proposed BEM technique based on CTM is faster and more accurate than FDM and FEM. Also, the proposed technique can be used in a wide variety of applications.
Nowadays, the knowledge of nonlinear thermal stresses in micropolar structures with temperature-dependent properties can be utilized by mechanical engineers in modeling of the solid composites, porous media and suspensions. As well as for chemists to detect the chemical reactions including adsorption and desorption.