Analysis of three-dimensional heat conduction in functionally graded materials by using a hybrid numerical method

International Journal of Heat and Mass Transfer 145 (2019) 118771

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Analysis of three-dimensional heat conduction in functionally graded materials by using a hybrid numerical method Wenzhen Qu a,b,⇑, Chia-Ming Fan c,d, Yaoming Zhang b a

School of Mathematics and Statistics, Qingdao University, Qingdao 266071, PR China Institute of Applied Mathematics, Shandong University of Technology, Zibo 255049, PR China c Department of Harbor and River Engineering & Computation and Simulation Center, National Taiwan Ocean University, Keelung 20224, Taiwan d Center of Excellence for Ocean Engineering, National Taiwan Ocean University, Keelung 20224, Taiwan b

a r t i c l e

i n f o

Article history: Received 24 May 2019 Received in revised form 11 August 2019 Accepted 21 September 2019

Keywords: 3D heat conduction Functionally graded materials Generalized finite difference method Krylov deferred correction method Long-time simulation

a b s t r a c t A hybrid numerical method is developed for three-dimensional (3D) heat conduction in functionally graded materials (FGMs) by integrating the advantages of the generalized finite difference method (GFDM) and Krylov deferred correction (KDC) technique. The temporal direction of the problems is first discretized by applying the KDC approach for high-accuracy results, which yields a partial differential equation (PDE) boundary value problem at each time step. The corresponding PDE boundary value problem is then solved by introducing the meshless GFDM that has no requirement of time-consuming meshing generation and numerical integration for 3D problems with complex geometries. Numerical experiments with four types of material gradations are provided to verify the developed combination scheme, and numerical results demonstrate that the method has a great potential for 3D transient heat conduction in FGMs especially for those in a long-time simulation. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction With the development of material processing technology, many material systems are widely designed and manufactured [1–4]. As a relatively new generation of composites, functionally graded materials (FGMs) are non-uniform microstructures with graded macro-properties (e.g. heat conductivity, specific heat, density, etc.) varied as a smooth function of spatial coordinates [5–7]. Due to the superior material properties, FGMs offer the benefit for resisting severe thermal loading. Thus, it is indispensable to analyze thermal properties of FGMs. In the past decades, many efforts have been devoted to investigations related with heat conduction analysis in FGMs by using the numerical approaches, such as the finite element method (FEM) [8–11], the finite difference method (FDM) [12], the boundary element method (BEM) [13– 16], and the meshless method [17–23]. Recently, many researchers pay attentions to the study of meshless methods which construct the trial functions in the scattered nodes [24–27]. We focus on strong-form meshless methods in this paper because of their being free from mesh generation and numerical quadrature [20–22,28,29]. Marin [28] applied the ⇑ Corresponding author at: School of Mathematics and Statistics, Qingdao University, Qingdao 266071, PR China. E-mail address: [email protected] (W. Qu). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118771 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

method of fundamental solutions (MFS) to the Cauchy problem for steady-state heat conduction in two-dimensional (2D) FGMs. Li et al. [21] extended the MFS for the solution of threedimensional (3D) transient heat conduction problems in FGMs with thermal conductivity and specific heat vary exponentially in z-direction. Skouras et al. [20] introduced the strong form meshless point collocation (MPC) method for heat conduction problems in isotropic and FGMs. Fu et al. [22] combined the Laplace transformation and boundary knot method (BKM) to investigate thermal behavior of slender FGMs with exponential variations. In above-mentioned strong-form meshless methods, the GFDM [30,31] is very simple, efficient, and easy to implement, which constructs explicit formulae for partial derivatives of undetermined variables by using Taylor series expansions and moving-least squares (MLS) approximations. Then enforcing the satisfactions of governing equation at interior nodes and boundary conditions at boundary nodes yield a system of linear or non-linear algebraic equations. Up to now, the GFDM has been applied to the solution of various engineering applications [32–36], especially for problems involving complicated geometry and high dimensions in analysis of non-linear and optimization problems [37–39]. In this paper, we present a numerical framework for the solution of 3D transient heat conduction in FGMs by introducing the GFDM with KDC technique. The time-dependent variable of the governing differential equation is first discretized by applying the

2

W. Qu et al. / International Journal of Heat and Mass Transfer 145 (2019) 118771

Nomenclature B c D/ H n p q Q S t u U  U x

Greek symbols boundary boundary with known temperature boundary with known flux d error q mass density / unknown function in boundary value problem x weighting function X domain

residual function specific heat partial derivative vector implicit function unit outward normal vector number of Gaussian nodes flux on boundary heat source density spectral integration matrix time temperature time derivative of temperature provisional solution of U vector of coordinates

C Cu Cq

Subscripts e exact solution n Numerical solution

KDC approach [39–41]that is of arbitrary order of accuracy and very stable. In the process of discretization, the KDC uses the spectral deferred correction technique as a preconditioner for the collocation formulations and applies the Newton-Krylov technique to the solution of resulting system of equations. At each time step, the corresponding PDE boundary value problem arising in temporal discretization is then solved by the strong-form meshless GFDM with no requirement of time-consuming meshing generation and numerical integration for 3D problems with complex geometries. To verify the present method, we finally provide four numerical examples with graded macro-properties respectively varied as four kinds of elementary functions of spatial coordinates. The structure of the paper is organized as follows. In Section 2, we present the governing differential equation of 3D transient heat conduction problems in FGMs. Section 3 shows the details of the combined scheme of the KDC and the GFDM. In Section 4, we provide four 3D numerical experiments including one case with the complicated geometry. Section 5 summarizes numerical results and makes some discussions. 2. Transient heat conduction problems in FGMs (3D) The governing differential equation for 3D general transient heat conduction in a bounded domain X with boundary C is given as

r  ðkðxÞruðx;t ÞÞ þ Q ðx; tÞ ¼ qðxÞcðxÞ

@uðx; tÞ ; x 2 X; 0 6 t < 1; @t

ð1Þ

where x ¼ ðx; y; zÞ, uðx; t Þ is the temperature function, Q ðx; t Þ is the heat source density, kðxÞ is the thermal conductivity, qðxÞ is the mass density, and cðxÞ is the specific heat. Two kinds of boundary conditions are given as 

uðx; t Þ ¼ uðx; tÞ; x 2 Cu ; 0 6 t < 1; qðx; tÞ ¼ kðxÞ

ð2Þ

@uðx; t Þ  ¼ qðx; t Þ; x 2 Cq ; 0 6 t < 1; @n



ð3Þ 

in which uðx; t Þ is known temperature on boundary Cu , q denotes known flux on boundary Cq , n is the unit outward normal vector at x, and C ¼ Cu [ Cq . The initial condition at a specific time t0 is given as

uðx; t 0 Þ ¼ u0 ðxÞ; x 2 X:

ð4Þ

3. The combined approach This section introduces two techniques of the KDC and the GFDM which serve as the fundamental building blocks for a combined algorithm used for 3D transient heat conduction problems in FGMs. We provide the combination details of the developed method. 3.1. Krylov deferred correction (KDC) ðx;t Þ We use U ðx; tÞ ¼ @u@t as the new unknown variable to remove the numerically unstable differential operator, which is similar to spectral deferred correction methods in [40]. Eq. (1) is then rewritten as

    Z t Z t kðxÞr2 u0 þ U ðx; sÞds þ rkðxÞr u0 þ U ðx; sÞds 0

0

þ Q ðx; t Þ ¼ qðxÞcðxÞU ðx; tÞ; x 2 X;

ð5Þ

where u0 ¼ uðx; 0Þ denotes the initial condition. To solve Eq. (5) as the computational time from t ¼ 0 to t ¼ T f , we usually separate   the time domain 0; T f into sub-intervals by using node points     T 0 ¼ 0; T 1 ; :::; T n ¼ T f . Each sub-interval T j ; T jþ1 is called as one ‘‘big” time step, and its step size is Dt ¼ T jþ1  T j . We show the details of marching Eq. (5) in ½T 0 ; T 1  based on Gaussian type nodes,   and then calculate the solutions in T j ; T jþ1 ; ð j P 1Þ with a similar way after getting the results at T j . In ½0; Dt, the Gaussian type nodes  T are expressed as t ¼ t 1 ; t 2 ; :::; t p , and the temporal discretized   solution vector is then written as U ¼ U 1 ðxÞ; U 2 ðxÞ; :::; U p ðxÞ with U i ðxÞ ¼ U ðx; t i Þ. As shown in [41], the spectral integration matrix S is introduced to approximate the integral operator R ti U ð x; s Þd s  D tSU, which maps the function values at Gaussian 0 nodes to its integral at the same node points using orthogonal polynomial approximation and Gaussian integration. It should be pointed out that S is dense and precomputed integration matrix. With the help of the spectral integration matrix S, we can discretize Eq. (5) as

kðxÞr2 ðu0 þ DtSU Þ þ rkðxÞrðu0 þ DtSU Þ þ Q ¼ qðxÞcðxÞU;

ð6Þ

in which u0 ¼ ½u0 ðx; 0Þ; u0 ðx; 0Þ; :::; u0 ðx; 0ÞT denotes the initial con  T represents the funcdition, and Q ¼ Q ðx; t 1 Þ; Q ðx; t 2 Þ; :::; Q x; t p tion values of heat source density at Gaussian nodes ti ði ¼ 1; 2; :::; pÞ.

3

W. Qu et al. / International Journal of Heat and Mass Transfer 145 (2019) 118771

To directly solve Eq. (6), it is computationally inefficient due to the dense matrix S particularly for that with a large number p of Gaussian nodes. Instead, we first calculate a provisional solution h iT    U ¼ U 1 ðxÞ; U 2 ðxÞ; :::; U p ðxÞ by using a low-order method to Eq. (5),

and then can define an error   T d ¼ U  U ¼ d1 ðxÞ; d2 ðxÞ; :::; dp ðxÞ . Assume that U ðx; tÞ and dðx; t Þ 



are the continuous counterparts of U and d which can be expressed with the polynomial interpolation at Gaussian nodes. Replacing 

U ðx; t Þ with the sum of U ðx; t Þ and dðx; tÞ, Eq. (5) can be recast as

i o Rt h kðxÞr u0 þ 0 U ðx; sÞ þ dðx; sÞ ds i o n R t h þrkðxÞr u0 þ 0 U ðx; sÞ þ dðx; sÞ ds þ Q ðx; tÞ h i ¼ qðxÞcðxÞ U ðx; t Þ þ dðx; t Þ ; x 2 X: 2

n



Since U is known,

R ti



0

ð7Þ



the error d both satisfy the following PDE boundary value problem at each Gaussian node ti

kðxÞr2 /ðxÞ þ rkðxÞr/ðxÞ 

qðxÞcðxÞ Dt 0

/ðxÞ ¼ f ðxÞ; x 2 X;

0

...

0

Dt 2 Dt 2 .. . Dt 2 Dt 2

0

...

0

Dt 3 .. . Dt 3 Dt 3

... .. .

0 .. .

. . . Dtp1 . . . Dtp1

0

3

0 7 7 7 0 7 7 ; .. 7 . 7 7 7 0 5

ð8Þ

Dt p

where Dt i ¼ ti  t i1 ði ¼ 1; :::; p; t 0 ¼ 0Þ represents the ‘‘small” time step size. Based on the matrix notation, we can evaluate the low 

order approximation d of the error d by the following formula:

          kðxÞr2 u0 þ DtS U þDt S d þ rkðxÞr u0 þ DtS U þDt S d þ Q

  ¼ qðxÞcðxÞ U þ d :



wðxÞ ¼ kðxÞ

ð13Þ

@/ðxÞ  ¼ wðxÞ; x 2 Cw ; @n

ð14Þ 





ð9Þ

  H U ¼ d;



Gaussian node ti ), Dt 0 ¼ t i  t i1 , f ðxÞ, /ðxÞ, and wðxÞ are all known 

at Gaussian node t i . It is worthy to note that: (i) U i ðxÞ and its normal gradient on the boundary can be respectively obtained based on Eq. 

(2) and Eq. (3); (ii) both di ðxÞ and its normal gradient on the boundary are set to be zero. That is to say, we do not directly use the temperature and its normal gradient but their time derivatives in the implementation of the developed method. The GFDM approach is applied to the solution of the boundary value problem of Eqs. (12)–(14). As shown in Fig. 1, we first distribute an irregular cloud of collocation points in the computational domain X. For each node x0 , treated as the central node, we can define a distance dm from the central node x0 , and call m nearest nodes xi ði ¼ 1; 2; :::; mÞ (i.e., jxi  x0 j 6 dm ) as the supporting nodes or neighbors. For each central node, its corresponding supporting nodes are defined as a ‘star’, which indicates that each node in X has its associated star. We then set /0 to be the value of /ðxÞ at the central node x0 and /i ði ¼ 1; 2; :::; mÞ to be function values at m nearest nodes xi ði ¼ 1; 2; :::; mÞ in the star. Based on the Taylor series expansion, the values of /i around the central node x0 can be expanded as [43,44]

 2 2 2 2 2 2 0 0 0 /i ¼ /0 þ hi @/ þ ki @/ þ li @/ þ 12 hi @@x/20 þ ki @@x/20 þ li @@x/20 @x1 @x2 @x3 1

2

In the KDC approach, Eq. (9) is considered as the following implicit function:

2

2

þhi ki @x@ 1/@x02 þ hi li @x@ 1 /@x03 þ ki li @x@ 2 /@x03 þ ::::::;

2

3

i ¼ 1; 2; :::; m: ð15Þ

where hi ; ki ; li can be respectively expressed as

ð10Þ 

where input and output variables are the provisional solution U and 



the error d respectively. The output variable d is equal to zero once 

the input variable U satisfies the collocation formulation in Eq. (6), i.e., Eq. (6) is equivalent to the following form:

 H U ¼ 0:

ð12Þ

with the boundary conditions





U ðx; sÞds can be approximated by using

0









Dt 1



approach, components U i of the provisional solution U and di of

where /ðxÞ represents U i ðxÞ or di ðxÞ (i.e., U ðx; t Þ or dðx; t Þ at each

based rectangular rule, this work gives S of the resulting backward Euler’s method as

Dt 1 6 Dt 6 1 6 6 Dt 1  6 Dt S ¼ 6 . 6 .. 6 6 4 Dt 1

After temporal discretization of the problems with the KDC

/ðxÞ ¼ /ðxÞ; x 2 C/ ;

DtS U with the high order spectral integration matrix S, i.e.,  R ti  U ðx; sÞds  DtS U . For the unknown error d, we use a low order 0   Rt quadrature to approximate 0i dðx; sÞds  Dt S d, and S represents a low order integration matrix. By applying the right-end point

2

3.2. Generalized finite difference method (GFDM)

ð11Þ

By applying the implicit function theorem to Eq. (9), we can

 observe that the Jacobian matrix J H of H U ¼ 0 is better conditioned than that of the previous collocation formulation in Eq. (6) because of the preconditioner as the result of low-order deferred correction procedure. The Jacobian-free Newton-Krylov (JFNK) technique is very suitable for solving such well-conditioned systems, thus the KDC method directly applies existing JFNK solvers [42] to the solution of the preconditioned system in Eq. (11).

Fig. 1. An irregular cloud of points and the selection of stars.

4

W. Qu et al. / International Journal of Heat and Mass Transfer 145 (2019) 118771

hi ¼ x1i  x10 ; ki ¼ x2i  x20 ; li ¼ x3i  x30 :

ð16Þ

  by denoting x0 ¼ x10 ; x20 ; x30 and xi ¼ x1i ; x2i ; x3i By truncating the Taylor series (15) after the second-order derivatives, a residual function Bð/Þ is defined as m h

P h2 0 0 0 /0  /i þ hi @/ þ ki @/ þ li @/ þ 2i @x1 @x2 @x3

Bð/Þ ¼

i¼1

@ 2 /0 @x21

þ

k2i @ 2 /0 2 @x2 2

l2 @ 2 / 0 @x23

ð17Þ in which xi ¼ xðhi ; ki ; li Þ represents the weighting function at point xi . The weighting function xi in this paper is set to be

xi ¼ 1  6

di dm

2

 3  4 di di þ8 3 ; dm dm

di 6 dm ;

" #T @/0 @/0 @/0 @ 2 /0 @ 2 /0 @ 2 /0 @ 2 /0 @ 2 /0 @ 2 /0 ; ; ; ; ; ; ; ; ; @x1 @x2 @x3 @x21 @x22 @x23 @x1 @x2 @x1 @x3 @x2 @x3

Minimizing the residual function Bð/Þ with respect to the partial derivatives D/ , i.e.,

ð20Þ

we then obtain the following linear equation system at each node

AD/ ¼ b;

ð21Þ

where

A¼

m X

a

 ðiÞ ð iÞ diag E1 diag E2 ;

ð22Þ

i¼1

is a symmetric matrix that stores geometric information for every ‘ node inside the star, represents a 9  9 matrix with each element



 ðiÞ ðiÞ equals to one, diag E1 and diag E2 are two diagonal matrices with diagonal elements respectively defined as

h ð iÞ E 1 ¼ xi h i ð iÞ

E2 ¼

xi h 2

2hi

ki

li

hi

2

ki

2

li

2ki

2li

hi

2

2

ki

2

iT hi k i 2

li

hi li

ki li

ð23Þ

; iT

2hi ki

2hi li

2ki li

;

ð24Þ

and b is the 9  1 vector expressed as 0

/0

/m

1

where E ¼ A B. Through the above-mentioned process, as shown in Eq. (26), the derivatives of the unknown function / at x0 are expressed as a linear combination of function values at its neighboring nodes. Then, the left side of governing Eq. (12) at node x0 can be recast as

kðx0 Þr2 /0 þ rkðx0 Þr/0 

qðx0 Þcðx0 Þ Dt

0

/0 ¼ m0 /0 þ

m X

mi /i :

ð27Þ

i¼1

The coefficients mi ði ¼ 0; 1; :::; mÞ in Eq. (27) can be easily determined by using Eq. (26). If x0 is a boundary collocation point with Dirichlet boundary condition, one can have 

/0 ¼ /ðx0 Þ:

ð28Þ

kðx0 Þ

m X @/0 nj ¼ m00 /0 þ m0 i /i ; xj i¼1

;

ðmþ1Þ1

where U ¼ ½/0 ; /1 ; /2 ; . . . ; /m  is a vector constructed by function values at all nodes inside the star. We refer interested readers to [43,44] for the details of the explicit expressions of A and b. With the help of Eqs. (21)–(25), the partial derivative vector D/ is recast as

j ¼ 1; 2; 3;

ð29Þ

where nj ð j ¼ 1; 2; 3Þ are components of normal vector n at x0 , coefficients m0i ði ¼ 0; 1; :::; mÞ can be easily determined by using the first-order partial derivatives in D/ of Eq. (26). If x0 is a boundary collocation point with Robin boundary condition, we can directly use the combined formulation of Eqs. (28) and (29). We can finally construct a system of linear algebraic equations by repeating the above procedures to each collocation point. Then we can get numerical results of the function / at collocation points after solving the system of equations. Fortunately, a sparse matrix system arises in the discretization of the GFDM, thus we can quickly and efficiently solve the final linear algebraic equations. 4. Numerical experiments We provide four numerical examples in this section. To fully verify the accuracy and convergency of the present method, graded macro-properties (the thermal conductivity, the mass density, and the specific heat) of materials are chosen as different elementary functions or their combinations in these examples. We test the accuracy of the numerical results by using the following two types of errors [45–48]:

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi XM  2 XM  2 ukn  uke = uke ; k¼1 k¼1

Max Error ¼ k un  ue k1

ð25Þ T

ð26Þ

/m

Global Error ¼

B C B C B /1 C B C B C ! m B C X B /2 C b ¼ BU ¼  xi Eð2iÞ ; x1 Eð21Þ ; x2 Eð22Þ ;  ; xm Eð2mÞ B C B C i¼1 C 9ðmþ1Þ B B . C B .. C B C @ A

C C C C C; C C A

If x0 is a boundary collocation point with Neumann boundary condition, we can obtain

ð19Þ

@Bð/Þ   ¼ 0; @ D/

1

1

ð18Þ

where di ¼ k xi  x0 k, and dm represents the distance from x0 to the farthest node inside the star. To express the PDE of Eq. (12) at the central node x0 , we define a partial derivative vector as follows:

D/ ¼

/0

B/ B 1 B B 1 1 D/ ¼ A b ¼ A BU ¼ EB /2 B . B . @ .

þ 2i

 i2 2 2 2 þhi ki @x@ 1/@x02 þ hi li @x@ 1/@x03 þ ki li @x@ 2/@x03 xi ;



0

ð30Þ ð31Þ

where M is the number of collocation points in the domain, uke ðk ¼ 1; 2; :::; M Þ and ukn ðk ¼ 1; 2; :::; MÞ are respectively exact and    numerical solutions, ue ¼ uke k ¼ 1; 2; :::; M and un ¼ ukn k ¼ 1; 2; :::; Mg. It should be noted that M is less than the total number of collocation points because of excluding boundary collocation points. For all numerical examples, the number of supporting nodes for the stars is set to be m ¼ 25, and the error tolerance of JFNK solver in the KDC approach is chosen to be 108 . The proposed method is currently implemented in Matlab. The codes of the present method were run on a desktop with Intel Core i7-6700 3.40 GHz Processor with 16 GB RAM. In addition, we use the software HyperMesh to get the coordinates of collocation points for numerical examples.

5

W. Qu et al. / International Journal of Heat and Mass Transfer 145 (2019) 118771

4.1. Test problem 1: Cube with exponential material gradation As the first example, we consider a unit cube with constant temperature on two sides. The material property varies only in z-direction. The mixed boundary condition is given as

8   > < uðx; y; 0; tÞ ¼ 0 C; uðx; y; 1; tÞ ¼ 100 C; 2 qðx; 0; z; t Þ ¼ 0 w=m ; qðx; 1; z; t Þ ¼ 0 w=m2 ; > : qð0; y; z; tÞ ¼ 0 w=m2 ; qð1; y; z; t Þ ¼ 0 w=m2 :

ð32Þ

kðx; y; zÞ ¼ 5e2z ; cðx; y; zÞ ¼ ez ; qðx; y; zÞ ¼ ez :

ð33Þ

 The initial condition is uðx; y; z; 0Þ ¼ 100sin p2z . The thermal conductivity kðx; y; zÞ, the specific heat cðx; y; zÞ and the mass density qðx; y; zÞ are given as

The heat source density is prescribed as

5p2 sint þ cost sinðpzÞ  10psintcosðpzÞ

pz

pzo þ125p2 sin  500pcos : 2 2

Q ðx; y; z; t Þ ¼ e2z



The analytical solution of temperature is

uðx; y; z; t Þ ¼ sintsinðpzÞ þ 100sin

Fig. 2. Distribution of temperature at plane x ¼ 0:5.

pz 2

ð35Þ

We run the simulation from t ¼ 0 s to t ¼ 10 s with a large time step size Dt ¼ 1 s. Fig. 2 shows the temperature distribution of the analytical solution at plane x ¼ 0:5 when time is t ¼ 10 s. The developed method employs four kinds of combinations for Gaussian node number p in each time step size and collocation point number n in space, which are (a)p ¼ 2; n ¼ 81; (b)p ¼ 3; n ¼ 637; (c)p ¼ 4; n ¼ 1215; (d) p ¼ 5; n ¼ 6647 respectively. Fig. 3 plots contours of relative errors of the temperature u

Fig. 3. Contours of relative errors of temperature at plane x ¼ 0:5.

6

W. Qu et al. / International Journal of Heat and Mass Transfer 145 (2019) 118771

The thermal conductivity and the specific heat are

kðx; y; zÞ ¼ x2 þ y2 þ z2 ; cðx; y; zÞ ¼ ðx þ y þ zÞ2

ð38Þ

The mass density is a constant q ¼ 1 kg=m . The heat source density is prescribed as n o  2 Q ðx;y;z;tÞ ¼ ecost exþyþz ðx þ y þ zÞ sint þ 3 x2 þ y2 þ z2 þ 2ðx þ y þ zÞ 3

ð39Þ

Finally, the analytical solution is determined as

uðx; y; z; tÞ ¼ ecost exþyþz Fig. 4. Variation of global errors of temperature at all collocation points as a function of time step number.

at t ¼ 10 s and collocation points (not including boundary collocation points) inside the plane x ¼ 0:5. It can be observed from the figure that the present method obtains the satisfied results even using 81 collocation points and 2 Gaussian nodes. As we can also see, the distribution of collocation points influences on the location of the maximum relative error, but the error at each collocation point is rapidly decreased with increasing numbers p and n. We also investigate the feature of the developed method for a long-time simulation from t ¼ 0 s to t ¼ 100 s with time step size Dt ¼ 1 s. 4 Gaussian nodes and 1215 collocation points are used in this simulation. Fig. 4 provides the variation of global errors of temperature at all collocation points in the domain as a function of time step number. As we can see, the global errors are basically stable when marching the solution from t ¼ 0 s to t ¼ 100 s. 4.2. Test problem 2: Beam structure with quadratic material gradation For the second example, we consider an FGM beam structure as shown in Fig. 5(a). The dimension of this structure is 0:50 m  0:6 m  0:4 m. The material property varies in three dimensions. The boundary and initial conditions are respectively assumed as

   

      u x; y; z; t ¼ ecost ex þ y þ z ; x; y; z 2 C

ð36Þ

uðx; y; z; 0Þ ¼ exþyþzþ1

ð37Þ

ð40Þ

For this case, we distribute 26,401 collocation points as shown in Fig. 5(b). The present method is used to march the solution from t ¼ 0 s to t ¼ 5 s using time step size Dt ¼ 1 s, and 5 Gaussian nodes are adopted in each time step size. Fig. 6 plots the relative @u ; @u ; @u , and @n on the surface of beam strucerrors of heat fluxes @u @x @y @z ture at t ¼ 5 s. It can be found that the accurate results of heat fluxes on the boundary are obtained by the proposed method. Next, we consider a long-time simulation for the present method with different number of Gaussian nodes in each time step size. Fig. 7 plots the temperature distribution of the analytical solution at plane z ¼ 0:4 when time is t ¼ 50 s. Fig. 8 shows the variations of global errors of temperature at all collocation points in the domain by using the developed method to march the solution from t ¼ 0 s to t ¼ 50 s with time step size Dt ¼ 1 s using 2, 3, and 4 Gaussian nodes. From the figure, we can see that the present method obtains lower relative errors with an increasing number of Gaussian points in each time step size. It also can be observed that these errors for different number of Gaussian nodes generally have no upward drift from 1 time step to 50 time steps, which is one of important features of the proposed method for long-time simulations. 4.3. Test problem 3: Serrated structure (thin body) with logarithmic material gradation For the third example, we consider the transient heat conduction in a serrated structure as shown in Fig. 9(a). The dimension of this structure is 1:0 m  0:3 m  0:033 m. The material property varies only in x-direction. The boundary and initial conditions are respectively given as

Fig. 5. The sketch of a beam structure (a) and the distribution of 26,401 collocation points (b).

W. Qu et al. / International Journal of Heat and Mass Transfer 145 (2019) 118771

7

Fig. 6. Relative errors of heat fluxes on the surface of beam structure.

Fig. 8. Variations of global errors of temperature at all collocation points in domain. Fig. 7. Distribution of temperature at plane z ¼ 0:4.

   

       u x; y; z; t ¼ cosð5tÞ x þey þ sin z ; x; y; z 2 C

The thermal conductivity is described as

ð41Þ

kðx; y; zÞ ¼ lnðx þ 2Þ

ð43Þ

The specific heat and the mass density are respectively assumed 

uðx; y; z; 0Þ ¼ x þ e þ sinz y

ð42Þ

to be constants c ¼ 2 J=kg C q ¼ 5 kg=m3 . The heat source density is

8

W. Qu et al. / International Journal of Heat and Mass Transfer 145 (2019) 118771

Fig. 9. The sketch of a serrated structure (a) and the distribution of 4961 collocation points (b).

4.4. Test problem 4: Electric motor (complicated domain) with trigonometric material gradation. The last numerical example is an electric motor as shown in Fig. 13(a), which has a dimension 0:40 m  1:00 m  0:63 m. The trigonometric material variations of the thermal conductivity kðx; y; zÞ, the specific heat cðx; y; zÞ and the mass density qðx; y; zÞ are respectively defined as

kðx; y; zÞ ¼ sinx þ cosy þ sinz; cðx; y; zÞ

Fig. 10. Distribution of temperature at plane z ¼ 0:033.

2

¼ sinðxyzÞ; qðx; y; zÞ ¼ sin z Q ðx; y; z; tÞ ¼ 50sinð5tÞðx þ e þ sinzÞ   1 þ lnðx þ 2Þðey  sinzÞ  cosð5t Þ xþ2 y

The boundary and initial conditions are respectively given as

ð44Þ

     3 3 3     u x; y; z; t ¼ sin et x þ y þ z ; x; y; z 2 C

ð47Þ

and

The analytical solution is finally given as

uðx; y; z; t Þ ¼ cosð5tÞðx þ ey þ sinzÞ

ð46Þ

 uðx; y; z; 0Þ ¼ sinð1Þ x3 þ y3 þ z3 ð45Þ

In the simulations, we adopt two kinds of collocation distributions: 4961 collocation points as shown in Fig. 9(b) and 14,852 collocation points. The temperature distribution of the analytical solution at plane z ¼ 0:033 and time t ¼ 5 s is given in Fig. 10. To march the solution from t ¼ 0 s to t ¼ 5 s, the present method uses time step size Dt ¼ 0:2 s and 5 Gaussian nodes in each time step @u on the sursize. Fig. 11 plots the relative errors of normal fluxes @n face of the serrated structure at t ¼ 5 s. We can find that the present method obtains accurate numerical results. Moreover, the relative errors of normal fluxes decay quickly when increasing number of collocation points from 4961 to 14,852. The stability of numerical errors is also investigated for the developed method in a long-time simulation from t ¼ 0 s to t ¼ 20 s with time step size Dt ¼ 0:2 s. 5 Gaussian nodes and 4961 collocation points are used in this simulation. Fig. 12 shows variations of global and max errors of temperature u at all interior collocation points with an increasing number of time step. It is obvious that the errors have no remarkable increase.

ð48Þ

The heat source density is

Fig. 12. Variations of global and max errors of temperature at all interior collocation points.

Fig. 11. Relative errors of normal fluxes on the surface of the serrated structure.

W. Qu et al. / International Journal of Heat and Mass Transfer 145 (2019) 118771

9

Fig. 13. The sketch of an electric motor (a) and the temperature distribution on the boundary (b).

Fig. 14. Global error variation of temperature at all interior collocation points with different time step size.

 2 Q ðx; y; z; tÞ ¼ et cosðet Þsin zsinðxyzÞ x3 þ y3 þ z3   3sinðet Þ x2 cosx  y2 siny þ z2 cosz þ 2ðx þ y þ zÞðsinx þ cosy þ sinzÞ ð49Þ

The analytical solution is determined as

  uðx; y; z; t Þ ¼ sin et x3 þ y3 þ z3

ð50Þ

We distribute 95,821 collocation points in space. The simulation is first run from t ¼ 0 s to t ¼ 2 s by using the developed method with 5 Gaussian nodes in each time step size. Fig. 13(b) shows the temperature distribution of the analytical solution on the boundary at t ¼ 2 s. Fig. 14 provides the global error variation of the temperature at all interior collocation points as a function of time step size Dt. We can see that the error rapidly decays with a decreasing time step size. It also can be observed that the convergence rate of the present method is approximate to 3.7 as indicated by the trend line. Under the above-mentioned assumption, we rerun the simulation with the time step size Dt ¼ 0:2 s. Fig. 15 @u on the boundary at plots the relative errors of normal fluxes @n t ¼ 2 s. We can observe that the max error of normal fluxes is less than 1.0E02.

Fig. 15. Relative errors of normal fluxes on the boundary.

technique for temporal discretization and the meshless GFDM approach for spatial discretization. To verify the present method, we provide numerical results of four experiments with different material gradations, and summarize some advantages: (a) The error curves of numerical results in Figs. 8 and 12 obtained by the present method have no upward drift with an increasing number of time step, which is an important feature for long-time simulations; (b) The present method allows the large time step size such as Dt ¼ 1 s for the problem with the analytical solution invloving the time term sinðtÞ; (c) The convergence rate of the present method is high and approximated to 3.7 even for the numerical example 4 with the complicated domain. We apply the proposed method for the solution of relatively simple problems in this paper. For the method, optimized implementation and further analyses are necessary to sufficiently explore its stability and accuracy. It also will be reported in the near future that the present method is extended for nonlinear transient heat conduction problems with thermal conductivity and specific heat varied in both temporal and spatial directions. Declaration of Competing Interest The authors declare that they have no conflicts of interest in the authorship or publication of this contribution.

5. Conclusion and generalization Acknowledgements In this paper, we develop a highly accurate numerical algorithm for the solution of 3D transient heat conduction problems in FGMs especially for cases in long-time simulations, based on the KDC

The work described in this paper was supported by the National Natural Science Foundation of China (Nos. 11802165), the Natural

10

W. Qu et al. / International Journal of Heat and Mass Transfer 145 (2019) 118771

Science Foundation of Shandong Province of China (Grant Nos. ZR2017BA003, ZR2017MA021), and the China Postdoctoral Science Foundation (Grant No. 2019M650158). References [1] A. Chattopadhyay, A. Pattamatta, Energy transport across submicron porous structures: a Lattice Boltzmann study, Int. J. Heat Mass Transf. 72 (2014) 479– 488. [2] K. Wang, L. Li, W. Xue, S. Zhou, Y. Lan, H. Zhang, Z. Sui, Electrodeposition synthesis of PANI/MnO2/graphene composite materials and its electrochemical performance, Int. J. Electrochem. Sci. 12 (2017) 8306–8314. [3] J.P. Pascon, Large deformation analysis of functionally graded viscohyperelastic materials, Comput. Struct. 206 (2018) 90–108. [4] H. Nazir, M. Batool, F.J. Bolivar Osorio, M. Isaza-Ruiz, X. Xu, K. Vignarooban, P. Phelan, Inamuddin, A.M. Kannan, Recent developments in phase change materials for energy storage applications: a review, Int. J. Heat Mass Transf. 129 (2019) 491–523. [5] A. Shaker, W. Abdelrahman, M. Tawfik, E. Sadek, Stochastic finite element analysis of the free vibration of functionally graded material plates, Comput. Mech. 41 (5) (2008) 707–714. [6] D. Boussaa, Optimization of temperature-dependent functionally graded material bodies, Comput. Methods Appl. Mech. Eng. 198 (37) (2009) 2827– 2838. [7] J.H. Tian, K. Jiang, Heat conduction investigation of the functionally graded materials plates with variable gradient parameters under exponential heat source load, Int. J. Heat Mass Transf. 122 (2018) 22–30. [8] Q.-H. Qin, Trefftz finite element method and its applications, Appl. Mech. Rev. 58 (5) (2005) 316–337. [9] Z.J. Fu, Q.H. Qin, W. Chen, Hybrid-Trefftz finite element method for heat conduction in nonlinear functionally graded materials, Eng. Comput. 28 (5) (2011) 578–599. [10] P. Liu, T. Yu, T.Q. Bui, C. Zhang, Y. Xu, C.W. Lim, Transient thermal shock fracture analysis of functionally graded piezoelectric materials by the extended finite element method, Int. J. Solids Struct. 51 (11) (2014) 2167– 2182. [11] V.N. Burlayenko, H. Altenbach, T. Sadowski, S.D. Dimitrova, A. Bhaskar, Modelling functionally graded materials in heat transfer and thermal stress analysis by means of graded finite elements, Appl. Math. Model. 45 (2017) 422–438. [12] B.-L. Wang, Z.-H. Tian, Application of finite element–finite difference method to the determination of transient temperature field in functionally graded materials, Finite Elem. Anal. Des. 41 (4) (2005) 335–349. [13] J. Sladek, V. Sladek, C. Zhang, Transient heat conduction analysis in functionally graded materials by the meshless local boundary integral equation method, Comput. Mater. Sci. 28 (3) (2003) 494–504. [14] A.I. Abreu, A. Canelas, W.J. Mansur, A CQM-based BEM for transient heat conduction problems in homogeneous materials and FGMs, Appl. Math. Model. 37 (3) (2013) 776–792. [15] Y. Sun, Indirect boundary integral equation method for the Cauchy problem of the Laplace equation, J. Sci. Comput. 71 (2) (2017) 469–498. [16] W. Yao, B. Yu, X. Gao, Q. Gao, A precise integration boundary element method for solving transient heat conduction problems, Int. J. Heat Mass Transf. 78 (2014) 883–891. [17] H. Wang, Q.-H. Qin, Y.-L. Kang, A meshless model for transient heat conduction in functionally graded materials, Comput. Mech. 38 (1) (2006) 51–60. [18] J. Sladek, V. Sladek, C. Hellmich, J. Eberhardsteiner, Heat conduction analysis of 3-D axisymmetric and anisotropic FGM bodies by meshless local PetrovGalerkin method, Comput. Mech. 39 (3) (2007) 323–333. [19] A. Khosravifard, M.R. Hematiyan, L. Marin, Nonlinear transient heat conduction analysis of functionally graded materials in the presence of heat sources using an improved meshless radial point interpolation method, Appl. Math. Model. 35 (9) (2011) 4157–4174. [20] E.D. Skouras, G.C. Bourantas, V.C. Loukopoulos, G.C. Nikiforidis, Truly meshless localized type techniques for the steady-state heat conduction problems for isotropic and functionally graded materials, Eng. Anal. Boundary Elem. 35 (3) (2011) 452–464. [21] M. Li, C.S. Chen, C.C. Chu, D.L. Young, Transient 3D heat conduction in functionally graded materials by the method of fundamental solutions, Eng. Anal. Boundary Elem. 45 (2014) 62–67. [22] Z.-J. Fu, Q. Xi, W. Chen, A.H.D. Cheng, A boundary-type meshless solver for transient heat conduction analysis of slender functionally graded materials with exponential variations, Comput. Math. Appl. 76 (4) (2018) 760–773.

[23] F. Wang, Q. Hua, C.-S. Liu, Boundary function method for inverse geometry problem in two-dimensional anisotropic heat conduction equation, Appl. Math. Lett. 84 (2018) 130–136. [24] X. Li, J. Zhu, On a Galerkin boundary node method for potential problems, Adv. Eng. Softw. 42 (11) (2011) 927–933. [25] W. Qu, W. Chen, Solution of two-dimensional stokes flow problems using improved singular boundary method, Adv. Appl. Math. Mech. 7 (1) (2015) 13– 30. [26] X. Li, Three-dimensional complex variable element-free Galerkin method, Appl. Math. Model. 63 (2018) 148–171. [27] Y. Gu, C.-M. Fan, R.-P. Xu, Localized method of fundamental solutions for largescale modeling of two-dimensional elasticity problems, Appl. Math. Lett. 93 (2019) 8–14. [28] L. Marin, Numerical solution of the Cauchy problem for steady-state heat transfer in two-dimensional functionally graded materials, Int. J. Solids Struct. 42 (15) (2005) 4338–4351. [29] J. Lin, C. Zhang, L. Sun, J. Lu, Simulation of seismic wave scattering by embedded cavities in an elastic half-plane using the novel singular boundary method, Adv. Appl. Math. Mech. 10 (2018) 322–342. [30] T. Liszka, An interpolation method for an irregular net of nodes, Int. J. Numer. Meth. Eng. 20 (9) (1984) 1599–1612. [31] J. Orkisz, P. Lezanski, P. Przybylski, Multigrid Approach to Adaptive Analysis of B. V. Problems by The Meshless GFDM, Springer, Netherlands, Dordrecht, 1999, pp. 173–180. [32] J.J. Benito, F. Ureña, L. Gavete, Influence of several factors in the generalized finite difference method, Appl. Math. Model. 25 (12) (2001) 1039–1053. [33] J.J. Benito, F. Ureña, L. Gavete, R. Alvarez, An h-adaptive method in the generalized finite differences, Comput. Methods Appl. Mech. Eng. 192 (5) (2003) 735–759. [34] P.-W. Li, W. Chen, Z.-J. Fu, C.-M. Fan, Generalized finite difference method for solving the double-diffusive natural convection in fluid-saturated porous media, Eng. Anal. Boundary Elem. 95 (2018) 175–186. [35] W. Qu, A high accuracy method for long-time evolution of acoustic wave equation, Appl. Math. Lett. 98 (2019) 135–141. [36] A. Shojaei, U. Galvanetto, T. Rabczuk, A. Jenabi, M. Zaccariotto, A generalized finite difference method based on the Peridynamic differential operator for the solution of problems in bounded and unbounded domains, Comput. Methods Appl. Mech. Eng. 343 (1) (2019) 100–126. [37] H.-F. Chan, C.-M. Fan, C.-W. Kuo, Generalized finite difference method for solving two-dimensional non-linear obstacle problems, Eng. Anal. Boundary Elem. 37 (9) (2013) 1189–1196. [38] S. Milewski, J. Orkisz, In search of optimal acceleration approach to iterative solution methods of simultaneous algebraic equations, Comput. Math. Appl. 68 (3) (2014) 101–117. [39] J. Chen, Y. Gu, M. Wang, W. Chen, L. Liu, Application of the generalized finite difference method to three-dimensional transient electromagnetic problems, Eng. Anal. Boundary Elem. 92 (2018) 257–266. [40] A. Dutt, L. Greengard, V. Rokhlin, Spectral deferred correction methods for ordinary differential equations, BIT Numer. Math. 40 (2) (2000) 241–266. [41] J. Huang, J. Jia, M. Minion, Arbitrary order Krylov deferred correction methods for differential algebraic equations, J. Comput. Phys. 221 (2) (2007) 739–760. [42] D.A. Knoll, D.E. Keyes, Jacobian-free Newton-Krylov methods: a survey of approaches and applications, J. Comput. Phys. 193 (2) (2004) 357–397. [43] L. Gavete, J.J. Benito, F. Ureña, Generalized finite differences for solving 3D elliptic and parabolic equations, Appl. Math. Model. 40 (2) (2016) 955–965. [44] Q. Hua, Y. Gu, W. Qu, W. Chen, C. Zhang, A meshless generalized finite difference method for inverse Cauchy problems associated with threedimensional inhomogeneous Helmholtz-type equations, Eng. Anal. Boundary Elem. 82 (2017) 162–171. [45] W. Qu, W. Chen, Y. Gu, Fast multipole accelerated singular boundary method for the 3D Helmholtz equation in low frequency regime, Comput. Math. Appl. 70 (4) (2015) 679–690. [46] J. Li, W. Chen, A modified singular boundary method for three-dimensional high frequency acoustic wave problems, Appl. Math. Model. 54 (2018) 189– 201. [47] W. Qu, C. Zheng, Y. Zhang, Y. Gu, F. Wang, A wideband fast multipole accelerated singular boundary method for three-dimensional acoustic problems, Comput. Struct. 206 (2018) 82–89. [48] W. Qu, C.-M. Fan, Y. Gu, F. Wang, Analysis of three-dimensional interior acoustic fields by using the localized method of fundamental solutions, Appl. Math. Model. 76 (2019) 122–132.