Efficient and accurate method for 2D periodic structures based on the physical features of the transient heat conduction

International Journal of Thermal Sciences 127 (2018) 213–231

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Efficient and accurate method for 2D periodic structures based on the physical features of the transient heat conduction

T

Q. Gao∗, H.C. Cui State Key Laboratory of Structural Analysis for Industrial Equipment, International Center for Computational Mechanics, Department of Engineering Mechanics, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116023, PR China

A R T I C L E I N F O

A B S T R A C T

Keywords: Periodic structure Transient heat conduction Matrix exponential Precise integration method Superposition principle

An efficient and accurate method is developed to solve the transient heat conduction problem in two-dimensional (2D) periodic structures. For a 2D periodic structure, according to the physical features of the transient heat conduction, the periodic property of the structure, and the physical meaning of the matrix exponential, it is demonstrated that the matrix exponential for a reasonable time step is a sparse matrix containing many identical elements. Next, based on the superposition principle of linear systems and the algebraic structure of the matrix exponential, computation of the response of original 2D periodic structure is transformed into computation of the responses of small-scale models with several unit cells. Finally, the precise integration method (PIM) is used to compute the temperature responses of the small-scale models. The proposed method not only inherits the accuracy and stability of the PIM but also achieves significantly improved computational time and storage requirements. A series of numerical examples demonstrate that the proposed method is more efficient than the Crank-Nicholson method and provides highly precise solutions even with a larger time step.

1. Introduction Periodic structures have attracted tremendous attention and substantial research in recent years on account of their many excellent thermal qualities and mechanical properties, which offer important economic benefits. These structures are composed of many identical structural components (unit cells) that are assembled end-to-end in one, two or three directions to form the complete system, yielding one-dimensional (1D), two-dimensional (2D), or three-dimensional (3D) periodic structures. Periodic structures have been applied to thermal protection systems, large satellite antennae and space trusses. Hence, it is very important to analyse the thermal behaviour of the periodic structure in many practical engineering applications. The multilayer composite material composed of different materials in a cohesive layered final product, which is a type of periodic structure, has found many applications in the food and consumer electronics packaging industries [1]. Thus, knowledge of heat fluxes and temperature distributions in laminated composites is very important. Early work on heat conduction in laminated composites can be attributed to the literature [2,3]. Horvay et al. [4,5] dealt with the thermal behaviour of periodic laminated composites where matrix and filler were stacked parallel to the heat propagation direction. A continuum mixture theory of the transient heat conduction in multilayered composites was

∗

proposed by Nayfeh [6]. Auriault [7] determined effective parameters for macroscopic descriptions for static or transient heat conduction in periodic laminated composites using a homogenization method. Matysiak et al. [8–10] solved the heat conduction problem in micro-periodic composites using a homogenized model with micro-local parameters. Noor et al. [11] developed a 2D computational model for nonlinear heat transfer in multilayered composite panels. Woźniak et al. [12] proposed a refined averaged model to solve the heat conduction problem in a rigid heat conductor with micro-periodic composite material, and using a similar method. Matysiak et al. [13] dealt with a mixed-in-time transient heat conduction problem in a periodically two-layered nondeformable half-space. Ignacazk et al. [14] employed the refined averaged heat conduction theory to obtain the formal solution to the 1D initial boundary value problem in a periodically layered plate consisting of many homogeneous isotropic layers. Bufler [15] dealt with the stationary heat conduction in a layered solid composed of planar, isotropic or transversely isotropic layers. Gottfried [16] employed a multi-scale approach based on a homogenization method of a periodic composite to analyse the aerothermal behaviour within transpirationcooled plates. Kulchytsky et al. [17–19] dealt with the heat conduction problem in semi-infinite periodically laminated layers by using the classical heat conduction description and the homogenized model with micro-local parameters. Abdelal et al. [20] investigated the transient

Corresponding author. E-mail addresses: [email protected] (Q. Gao), [email protected] (H.C. Cui).

https://doi.org/10.1016/j.ijthermalsci.2018.01.006 Received 13 August 2017; Received in revised form 4 January 2018; Accepted 4 January 2018 1290-0729/ © 2018 Elsevier Masson SAS. All rights reserved.

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

Nomenclature

t Ω Γ Γ2 q h ρ T0 T cp k Q n I B (kx , ky) η m tk tCN tFPIM tR x , y, z x Γ1

Γ3 Text a ∇ T0 T C K P H Ψα b(kx , ky) A rα Tk TREF TFPIM eR

time [s] problem domain boundary of Ω Neumann boundary heat flux [W m−2 ] heat transfer coefficient [W m−2 K−1] density [kg m−3] initial temperature [K or °C ] temperature [K or °C] specific heat [J kg−1 K−1] thermal conductivity [W m−1 K−1] heat source density [W m−2 ] unit outward normal vector unit matrix the basic finite element model time step length [s] number of unit cells affected by excitations within one time step time kη [s] CPU times of the Crank-Nicholson method [s] CPU times of the proposed method [s] relative time Cartesian coordinates [m] space-variable vector Dirichlet boundary

Robin boundary temperature of ambient medium [K or °C] thermal diffusivity coefficient [m2 s−1] temperature gradient operator initial temperature vector nodal temperature vector heat capacity matrix heat conductivity matrix thermal force vector system coefficient matrix matrices related to external excitations the small-scale model matrix exponential external excitation vectors temperature vector at time tk the reference solution [K or °C] solution obtained using the proposed method [K or °C] relative error

Abbreviations PIM FEM PCG 1D 2D 3D

precise integration method finite element method preconditioning conjugate gradient one dimensional two dimensional three dimensional

solids that can account for the scale effect of unit cells by applying the mathematical homogenization method and developed a multi-scale multi-physics analysis method for simulating the thermo-mechanical coupled behaviour characterized by the proposed model. By using a new second-order two-scale method, Ma et al. [33] solved a static heat conductive problem in a periodic porous domain with radiation boundary conditions on the surface of cavities. Yang et al. [34–38] studied the heat transfer performance of periodic porous materials using a novel multiscale asymptotic expansion method. In addition to periodic layered composite materials and periodic porous materials, Dasgupta et al. [39] used the multiscale asymptotic analysis method to predict the orthotropic thermal conductivity of plain-weave fabric-reinforced composite laminates. Alzina [40] estimated overall thermal conductivity in a triaxially braided composite by using multiscale modeling at cryogenic temperatures. Ozdemir et al. [41] used the computational homogenization method to simulate heat transfer in heterogeneous materials with detailed microstructural constituents. Kouznetsova et al. [42] proposed a novel second-order computational homogenization approach for the multi-scale analysis of multi-phase materials. Based on the extended finite element method (XFEM) and the level set method (LSM), Wang et al. [43] proposed an approach for evaluating the effective thermal conductivity of both particle- and fibre-reinforced composites that contain arbitrarily shaped interfaces and exhibit both strong and weak discontinuities. It is of great importance to analyse the transient heat conduction of periodic structures in practical engineering areas. However, analytical solutions for periodic structures with complex geometry and boundary conditions are usually not available [44]. As seen from the above literature review, the research studies regarding the solution of transient heat conduction problems of periodic structures mainly focus on homogenization methods and multiscale methods. Based on averaging techniques, homogenization methods can give a crude approximation of the overall behaviour of the periodic structure. However, the homogenized solution is not sufficient to describe local temperature fluctuations. Therefore, types of multiscale methods for heat conduction problems in periodic structures have been proposed. Multiscale

heat conduction of rectangular variable stiffness laminates using the finite element method (FEM). Muliana et al. [21] formulated a twoscale homogenization scheme to predict effective thermal conductivities and field variables of fibre-reinforced plastic laminated composites subject to external thermal stimuli. Chung et al. [22] described two multiscale approaches employing the asymptotic analysis homogenization method to study nonlinear thermal heat conduction in composite materials with temperature-dependent conductivity. Matine et al. [23,24] studied heat conduction problems within a periodic layered composite using the multiscale asymptotic expansion method, and additional asymptotic expansion terms were introduced to correct the edge effects in the heat conduction. Wang et al. [25] proposed an efficient parallel multiscale numerical algorithm for a parabolic equation with rapidly oscillating coefficients representing heat conduction in periodic layered composite materials. Monteiro et al. [26] proposed a method combining a multi-scale approach and a model reduction technique based on proper orthogonal decomposition to solve highly nonlinear conduction problems in structures made up of periodic heterogeneous materials. Another type of periodic structure is porous material, which has many elegant qualities, such as low relative density, heat insulation, high heat resistance, and so on, and has been widely used in a variety of industrial products. Porous material tends to have many very small pores, which may be open and connected or closed. However, because of the regular or sufficiently regular heterogeneity, the porous material is often regarded as the periodic structure in the study. Considering the periodic porous material, some interesting work has been proposed in the past few years. Liu et al. [27] presented a homogenization-based multi-scale method for predicting the effective thermal conductivity of porous materials with radiation. Yusuke et al. developed a homogenization method for the thermal analysis of a porous medium with ellipsoidal pores [28], resin composites with ellipsoidal fillers [29] and packed beds [30]. Zhai et al. [31] discussed the initial boundary value problem for heat transfer in porous (lattice-type) materials using the homogenization method and the multiscale asymptotic method. Terada et al. [32] derived a two-scale thermo-mechanical model for porous 214

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

structures can also be solved by the standard FEM, which involves discretizing the spatial domain using a mesh to obtain the ordinary differential equation in the time domain and then using direct timestepping methods to solve the ordinary differential equation. However, to obtain a high-accuracy solution for a periodic structure with complex geometry and boundary conditions, the huge numbers of degrees of freedom are also required to capture the micro-scale behaviour, which leads directly to significant computational resources in terms of CPU power and memory.

asymptotic analysis methods have been used to solve the heat conduction problems of composites with periodic microstructures based on an asymptotic expansion of the unknown temperature field [45,46]. By considering higher-order terms in the expansion of the temperature field, the influence of microstructures on the macroscopic approximation can be captured more accurately. However, these methods can yield huge numbers of degrees of freedom and high computational cost depending on the macro- and micro-structural details [20,37]. As we already know, the transient heat conduction problems in periodic

Fig. 1. (a) The two-dimensional periodic structure; (b) the unit cell; and (c) the FEM mesh of the unit cell.

215

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

conduction and the periodic property of the structure, which acts to significantly reduce the computational cost and storage requirements even for large-scale problems. As mentioned previously, FEM has its drawbacks. Nevertheless, as compared to other numerical approaches, FEM still has numerous advantages [64,65]. This method is completely general with respect to arbitrary geometric shape and material properties, and the matrices obtained by FEM are symmetric and have narrow band-width property, which is useful for this research. Thus, in this paper, to obtain highaccuracy solutions more efficiently for transient heat conduction problems in 2D periodic structures, the FEM is adopted to discretize the problem domain, and a computational approach is proposed to improve the computational efficiency based on the precise integration method (PIM), the superposition principle of linear systems, the physical feature of transient heat conduction and the periodic property of the structure. The remainder of this paper is outlined as follows. Next, the mathematical model and the finite element model for the transient heat conduction problem in 2D periodic structures are presented. We will briefly discuss the PIM in Section 3. The algebraic structure of the matrix exponential corresponding to the 2D periodic structures is analysed in detail in Section 4. In Section 5, computation of the response of the original 2D periodic structure is transformed into computation of the responses of the small-scale models with a few unit cells,

Various approaches aimed at reducing computational effort have been proposed. Based on preconditioning techniques and speeding up the calculation of the matrix-vector product, iterative solvers are usually used to reduce the computational cost [47–52]. Serra et al. [53] proposed an iterative method to reduce the processing time of simulating 3D heat diffusion. Based on an element-by-element strategy, iterative methods with preconditioning are employed for nonlinear transient heat conduction analysis by Winget et al. [54]. In the past few decades, for solving the transient heat conduction problems efficiently, the proper orthogonal decomposition technique has been used to construct some reduced models with the finite difference method [55], the finite element method [56,57], the finite volume method [58], and the classical element free Galerkin method [59]. Another technique for reducing the computational cost of the element free Galerkin method is mass lumping procedure, which has also been developed for transient heat conduction problems successfully [60]. In addition, Kevin et al. [61] developed a parallel domain decomposition Laplace transform boundary element method algorithm for the solution of large-scale transient heat conduction problems. Feng et al. [62] presented a radial integration boundary element method for solving three-dimensional transient heat conduction problems, and Hitti et al. [63] used the zone method to reduce the CPU time greatly. The following paper will formulate a new approach with the physical feature of transient heat

Fig. 2. The finite element model corresponding to the 2D periodic structure.

216

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

a = (k y − 1)[Nx (dtb + dlr + di + dc ) + dlr + dc ] + (k x − 1)(dtb + dc ) + 1 b = a + dtb + 2dc − 1 c = a + (Nx − k x + 1)(dtb − dc ) + (k x − 1)(dlr − 2dc + di ) + dc d = c + 2dlr + di − 1 e = k y Nx (dlr + dc + dtb + di ) + k y (dlr + dc ) + (k x − 1)(dtb + dc ) + 1

which can reduce the computational efforts significantly. The detailed algorithm for determining the reasonable time step and the scale of the small-scale model is described in Section 6. The numerical examples are presented in Section 7, which demonstrate the efficiency and accuracy of the proposed method. Finally, the conclusions are given in Section 8.

f = e + (dtb + 2dc ) − 1 (6)

2. FEM for transient heat conduction in a 2D periodic structure

Note that, based on the above discrete strategy, once the heat capacity and heat conductivity matrices of a unit cell are obtained, according to the correspondence relationship between nodes of all unit cells, the heat capacity matrix and the heat conductivity matrix of the entire periodic structure can be obtained by the accumulation of the matrices of the unit cell and the modification according to the Robin boundary condition. Based on FEM theory, the governing equation of transient heat conduction can be written in the following form [65,66]

In this section, we will introduce the governing equation and the finite element model of transient heat conduction for the 2D periodic structure. The 2D periodic structure, as shown in Fig. 1(a), is formed by Nx × Ny unit cells (see Fig. 1(b)) in two directions, where Nx and Ny represent the numbers of unit cells in the horizontal and vertical directions, respectively. The governing equation for transient heat conduction in the 2D periodic structure can be written as

ρ (x) cp (x)

∂T (x, t ) − ∇ [k (x) ∇T (x, t )] = Q (x, t ) in ∂t

Ω × (0, t f )

CT˙ + KT = P, T (0) = T0

(1)

where k (x) , ρ (x) and cp (x) are, respectively, the thermal conductivity, mass density and specific heat of the medium, which are periodic functions in the two directions. T (x, t ) is the temperature at location x at time t ; (0, t f ) is the time interval; the function Q (x, t ) is the heat source density; Ω stands for the considered domain enclosed by the boundary Γ ; and Γ = Γ1 ∪ Γ2 ∪ Γ3 . On the boundary, the Dirichlet, Neumann and Robin boundary conditions are prescribed as follows

T=T

on Γ1 × (0, t f )

(2)

n⋅k (x) ∇T (x, t ) = q

on Γ2 × (0, t f )

(3)

n⋅k (x) ∇T (x, t ) = h (T − Text )

on Γ3 × (0, t f )

(4)

where T and q are the given temperature and heat flux on the boundary; the variables h and Text are the heat transfer coefficient and the temperature of surroundings, respectively; and n is the unit outward normal vector to the boundary. The boundary conditions of the 2D periodic structure are shown in Fig. 1(a), where the blue dotted lines, black solid lines and red dashed lines represent the Dirichlet, Neumann and Robin boundary conditions, respectively. The initial condition is

T (x, 0) = T0 (x)

in Ω

(7)

where T is the temperature vector, T˙ is the derivative of T with respect to time, T0 is the initial temperature vector, K is the heat conductivity matrix, C is the heat capacity matrix, and P is the thermal force vector. It should be noted that Eq. (7) denotes the governing equation that addresses the Dirichlet boundary condition in this paper. The subject of this research is linear transient heat conduction, and we assume that the heat capacity matrix C and the heat conductivity matrix K are not dependent on temperature, which is a reasonable assumption for actual engineering problems when changes in temperature are sufficiently small during the process. And it is well known that the matrices C and K are also only related to the physical parameters, geometrical shape and type and location of the boundary conditions, and are independent of the values of the given temperature T , the given heat flux q , the temperature of surroundings Text , and the heat source density Q (x, t ) . Alternately, T , q , Text and Q (x, t ) are converted into the thermal force vector P . Therefore, for FEM, once the thermal force vector P is obtained, the boundary conditions can be set to be homogenous, i.e., the finite element model shown in Fig. 2 gives the same governing equation (7). The above obvious discussion is prepared for the convenience of the analysis in Section 5. 3. PIM for linear transient heat conduction The most popular methods used to solve the ordinary differential equation (7) are direct time stepping methods. In particular, the θ -difference method [67–70] plays a dominant role in solving transient heat conduction problems. However, for these methods, a sufficiently small time step is required to achieve high precision or good stability. Therefore, significant computational resources in terms of CPU power and memory are usually required for large-scale problems. Another method used to solve ordinary differential equations is PIM. Because this method can avoid significant truncation and round-off errors, it has unconditional stability and can give high-precision numerical results at the integration points. Therefore, PIM has been widely applied to many fields, such as heat conduction [71–73], wave propagation [74–76], random vibration [77,78], and optimum control [79,80]. However, PIM can be quite time-consuming for large-scale structures. The purpose of this paper is to develop an efficient and accurate method based on PIM. Thus, PIM will be illustrated briefly in this section. Eq. (7) can be rewritten as

(5)

in which T0 (x) is the initial temperature of the 2D periodic structure. To make full use of the periodic property of the structure, the following strategy is presented to discretize the solution domain. Because the geometry and physical parameters for all unit cells are the same, the same finite element discretization can be used for all unit cells of the 2D periodic structure. To achieve this purpose, it is required that there be a one-to-one correspondence between the nodes located on the top (left) and bottom (right) edges of the finite element mesh of each unit cell, i.e., when the top (left) edge of a unit cell and its bottom (right) edge coincide together after translation, the positions of the nodes located on the top (left) edge are the same as those of the nodes located on the bottom (right) edge one-by-one. Assume that each unit cell has di internal nodes, dc nodes at the four corners, dtb + 2dc nodes on its top or bottom edge, and dlr + 2dc nodes on its left or right edge (see Fig. 1(c)); the labels of unit cells composing the 2D periodic structure are (1,1), (1,2), …, (Nx , Ny ) (see Fig. 2). Then, the nodes of unit cell (k x , k y ) can be numbered according to the following rules, i.e., the nodes located on the bottom edge are numbered from a to b , the nodes located on the top edge are numbered from e to f , and the other nodes are numbered from c to d , in which

T˙ = HT + f

(8)

in which

H = −C−1K, f = C−1P

(9)

For numerical integration, the time domain is divided into a series of time intervals of equal length η , i.e., 217

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

t0 = 0, t1 = η , …, tk = kη , …

(10)

heat conduction problems in 2D periodic structures. In the next section, we will illustrate the algebraic structure of the matrix exponential A and the matrices Ψα corresponding to the 2D periodic structure.

Let Tk = T(tk ) ; then, the solution of Eq. (8) at discrete times can be derived as

Tk + 1 = ATk +

∫0

η

exp[H (η − ξ )] f (tk + ξ )dξ

(11)

4. Algebraic structure of the matrix exponential corresponding to 2D periodic structures

(12)

In this section, the special algebraic structure of the matrix exponential will be analysed via the physical meaning of the matrix exponential, the physical feature of transient heat conduction, and the periodic property of the structure. The physical meaning of the matrix exponential is given firstly. If the thermal force vector P does not exist in Eq. (7), i.e., P = 0 , Eqs. (11) and (19) can be simplified as

in which A is the exponential of the matrix Hη, i.e.,

A = exp(Hη)

When Eq. (11) is used to compute transient heat conduction problems, the matrix exponential A and the integration in Eq. (11) need to be computed efficiently and accurately. PIM is based on an addition theorem and the idea of computing and storing the incremental part of the matrix exponential. The main process for computing the matrix exponential is as follows. For the given matrix H and a time step η, let

H′ = Hη /2N

Tk + 1 = ATk

The following physical meaning of the matrix exponential can be clearly shown from Eq. (22). Suppose that the initial temperature of the j−th node is unity, and that the initial temperatures of the other nodes are zero. Then, the temperature response of the i−th node at t = η is A(i, j ) . The physical feature of transient heat conduction can be illustrated as follows. For a homogeneous unbounded medium, when an instantaneous heat source is applied at point (x 0 , y0 , z 0) at time t = t0 , the temperature response at (x , y, z ) is given by the following expression [85,86]

(13)

in which N is an integer. By using the addition theorem of exponential function (2N algorithm), Eq. (12) can be rewritten as

A = [exp(H′)]2

N

(14)

η/2N

is very small, and thus the exponential of If N is large enough, the matrix H′ can be computed approximately by a Taylor series with q terms, i.e.,

exp(H′) ≈ I + H′ + (H′)2 /2! + ⋯+(H′)q / q!

(15)

To avoid the round-off error, exp(H′) is further divided into two parts, i.e.,

exp(H′) ≈ I + R 0,

R 0 = H′ + (H′)2 /2! + ⋯+(H′)q / q!

L (t , x , y, z ) =

(16)

(17)

(18)

When the matrix exponential is obtained, the next problem is how to compute the integration in Eq. (11). For the most used heat loads, the integration in Eq. (11) can be performed analytically. For example, if the heat load is a combination of linear terms and harmonic terms in a time step, i.e., P = r0 + r1t + r2 sin(ωt ) + r3cos (ωt ) , then Eq. (11) can be rewritten as 3

Tk + 1 = ATk +

∑ Ψα rα α=0

(19)

in which

Ψ0 = (I − A) K−1, Ψ1 = tk Ψ0 + E,

Ψ2 = X sin(ωtk ) + Y cos(ωtk ) Ψ3 = X cos(ωtk ) − Y sin(ωtk )

1 (4πaη)3/2

⎨ ⎩ 0,

exp ⎡− ⎣

(x − x 0)2 + (y − y0 )2 + (z − z 0)2 ⎤, 4aη

⎦

t > t0 t ≤ t0

where a = k /(ρc ) is the thermal diffusivity coefficient and η = t − t0 is the time interval. L (t , x , y, z ) is called the fundamental solution of the three-dimensional transient heat conduction equation. From Eq. (23), we can see that the temperature response at time t follows an exponential decay with the distance between the location (x , y, z ) and the point source. Therefore, the temperature response at a location that is further away from the point source is closer to zero. As far as a nonhomogeneous material is concerned, although to obtain the analytic fundamental solution is very difficult, it is clear that the distribution characteristic of its temperature response is similar to that of the homogeneous material. In addition, for a discrete system, as the number of elements is increased, the solution obtained from the FEM can gradually converge to the exact solution of the continuous system [64]. So the physical feature of its solution is similar to that in the continuous system. The sparsity of the matrix exponential is proved as follows. For convenience, let d kjx , ky denote the j−th node of the unit cell (k x , k y ) . According to the physical feature of transient heat conduction, if only a unit initial temperature is applied at the node d kjx , ky of the 2D periodic structure, after a time step η, the non-zero temperature response only exist near the node d kjx , ky , and those of the nodes that are far enough

Finally, the matrix exponential can be given by

A = I + RN

⎧

(23)

Then, the addition theorem is applied to the incremental part, i.e.,

for i = 1: N R i = 2R i − 1 + R i − 1 × R i − 1 end

(22)

(20)

and

away from the node d kjx , ky are zero. Based on the physical meaning of

E = [(I − A) H−1 + ηI] K−1 X = (ω2I + H2)−1 [ω sin(ωη) I + [A − cos(ωη) I] H] C−1 Y = (ω2I + H2)−1 [ω [A − cos(ωη) I] − sin(ωη) H] C−1

the matrix exponential, the elements of the d kjx , ky −th column in the matrix exponential A have the following features: the elements corresponding to the nodes that are near to the node d kjx , ky are significantly greater than zero, while the elements corresponding to the nodes that are far away from the node d kjx , ky are zero. Therefore, it can be concluded that the matrix exponential corresponding to a small time step η is a sparse matrix. To prove the repeatability of elements in the matrix exponential, based on the physical feature of transient heat conduction and the periodic property of the structure, we assume that the initial temperature can only affect m unit cells in any direction within one time step. Thus, when k x , qx ∈ (m , Nx − m) and k y, qy ∈ (m , Ny − m) , the

(21)

where I is the unit matrix. PIM has many numerical advantages, such as unconditional stability, high precision and a zero-amplitude rate of decay [81–84]. Additionally, the matrix exponential computed by PIM can be regarded as the exact solution on the computer. However, the computational effort of PIM is O (n3) , so it is a quite time-consuming task when the matrix H is large. Therefore, in this paper, a method that can take full advantage of the PIM and avoid its disadvantage is proposed to solve the transient 218

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

temperature responses caused by a unit initial temperature that is applied at nodes d kjx , ky or d qjx , qy are not affected by the boundary conditions of the 2D periodic structure (see Fig. 3). Note that the position of the node d kjx , ky located on the unit cell (k x , k y ) is the same as that of the node

same algebraic structure. In the next section, an efficient and accurate algorithm that can reduce the storage of zero elements and the computation of repetitive elements is proposed for solving transient heat conduction problems of 2D periodic structures.

d qjx , qy located on the unit cell (qx , qy ) , and the geometric and physical properties of all unit cells are identical. Thus, the non-zero temperature responses for both cases are the same by an appropriate translational motion. Based on the above analysis, the response of node d pi x , py (px ∈ [k x − m , k x + m]; py ∈ [k y − m , k y + m]) caused by the unit in-

5. Efficient and accurate method for transient heat conduction in 2D periodic structures In this section, an efficient and accurate method is presented based on PIM, the superposition principle in linear systems, the physical feature of transient heat conduction, and the periodic property of the structure. For convenience of understanding and analysis, we will demonstrate the main idea of the proposed method by a 2D periodic structure with Nx = 5 and Ny = 5 in subsection 5.1. Furthermore, according to the physical feature of transient heat conduction discussed in Section 4, the temperature response caused by the excitations in a unit cell is localized. Thus, for this example, we assume that the excitations can only affect one unit cell in any direction within one time step η, i.e., m = 1. In subsection 5.2, the same idea of the illustration example is generalized directly to the general 2D periodic structure.

itial excitation of node d kjx , ky is the same as the response of node

d pi x + qx − kx , py + qy − ky caused by the unit initial excitation of node d qjx , qy . Thus, according to the physical meaning of the matrix exponential, if k x , qx ∈ (m , Nx − m) and k y, qy ∈ (m , Ny − m) , then

(

)

(

A d pi x , py , d kjx , ky = A d pi x + qx − kx , py + qy − ky , d qjx , qy

)

(24)

Eq. (24) shows that there are many repetitive non-zero elements in the matrix exponential corresponding to 2D periodic structures. In this section, we can conclude that, for a reasonable time step, the matrix exponential corresponding to the 2D periodic structure is a sparse matrix with large numbers of repeatable elements. Based on a similar analysis, it can be easily proven that the matrices Ψα have the

Fig. 3. Schematic description for the repeatability of elements in the matrix exponential corresponding to the 2D periodic structure.

219

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

5.1. Main ideas of the proposed method According to the analysis given in the last paragraph of Section 2, the finite element model corresponding to a 2D periodic structure with 5 × 5 unit cells is shown in Fig. 4(a), where the boundary conditions are homogenous and the shadow domains represent the external excitation rα and the initial temperature Tk . Based on the superposition principle, the external excitation rα can be subdivided into 25 basic external excitations (k , k ) r α x y (k x , k y = 1, 2, ⋯, 5) according to the following principles: (1) the (k , k )

basic external excitation r α x y indicates that only the nodes located on the mesh of unit cell (k x , k y ) are applied by the external excitations and that there are no external excitations on the remaining unit cells; and (2) the original external excitation rα can be obtained by adding to(k , k ) gether all 25 basic external excitations r α x y , i.e., 5

rα =

5

∑ ∑

(k , ky )

rα x

(α = 0,1,2,3) (25)

kx = 1 ky = 1

In the same manner, the initial temperature Tk can be expressed as the (k , k ) sum of 25 basic initial temperatures T k x y (k x , k y = 1, 2, ⋯, 5), i.e., 5

Tk =

5 (k , ky )

Tk x

∑ ∑

(26)

kx = 1 ky = 1

Substituting Eqs. (25) and (26) into Eq. (19), we have 5

Tk + 1 =

5

3

∑ ∑ ⎛AT (kkx,ky) + ∑ Ψα r (αkx,ky) ⎞ α=0 kx = 1 ky = 1 ⎝ ⎠

Fig. 5. The relationship between the basic finite element model B (1,1) and the small-scale finite element model b(1,1) .

(27)

Based on the above analysis, the solution of the original finite element model shown in Fig. 4(a) can be obtained by the superposition of the solutions of the 25 basic finite element models B (kx , ky) with basic ex(k , k ) (k , k ) ternal excitations r α x y and T k x y (see Fig. 4(b)). Because we have assumed that the excitations only affect one unit cell in any direction, after a time step, the non-zero temperature responses of each basic finite element model B (kx , ky) caused by the basic (k , k ) (k , k ) excitations r α x y and T k x y only exist in the local domain which is composed of 3 × 3 unit cells around the unit cell (k x , k y ) . Based on this fact, the basic finite element model B (1,1) is considered, and its solution domain can be subdivided into two domains Ω1 and Ω2 , in which the domain Ω1 is composed of the unit cells (1,1) , (1,2) , (1,3) , (2,1) , (2,2) , (2,3) , (3,1) , (3,2) and (3,3) , and the domain Ω2 is composed of the remaining unit cells, as shown in Fig. 5. Since the external excitation and initial temperature are only applied on the unit cell (1,1) , the non-zero temperature responses only exist on the unit cells (1,1) , (1,2) , (2,1) and (2,2) (see Fig. 5). Thus, there is no heat flux exchanged between the domains Ω1 and Ω2 . The basic finite element B (1,1) can be decoupled into

two structures with domains Ω1 and Ω2 , respectively, and the arbitrary homogenous boundary conditions can be applied on the top and right boundaries of the domain Ω1, while the boundary conditions of the bottom and left boundaries of the domain Ω1 remain the same as those of the basic finite element model B (1,1) . Because the temperature responses for the structure with domain Ω2 are always zero within a time step, if the temperature response of the small-scale model with domain Ω1 is obtained, the temperature response of the basic finite element model B (1,1) can be obtained by the correspondence relationship between the nodes. Based on the same analysis, computation of the temperature responses of the 25 basic finite element models can be transformed into computation of the temperature responses of the 25 small-scale models. The schematic diagrams of the 25 small-scale models are shown in Fig. 6, where the blue dash-dotted lines represent the arbitrary homogenous boundary conditions and the red solid lines indicate the original homogenous boundary conditions of the basic finite element models. Fig. 4. The superposition principle for transient heat conduction in the 2D periodic structure with 5 × 5 unit cells: (a) the finite element model corresponding to the 2D periodic structure; and (b) the corresponding basic finite element models. The shadow domain represents the excitations, ∑ represents summation.

220

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

Fig. 6. The small-scale models corresponding to the 2D periodic structure with 5 × 5 unit cells.

For convenience, the small-scale model corresponding to the basic finite element model B (kx , ky) is denoted by b(kx , ky) . It is important to note that, for the nine small-scale models enclosed by the solid lines in Fig. 6, because arbitrary homogenous boundary conditions can be applied on the boundary of each model, we can apply the same homogenous boundary conditions for all these models. Thus, according to the periodic property of the structure, the nine small-scale

models are the same except that their excitations may be different. Thus, computation of the responses of the nine small-scale models can be transformed into computation of the responses of one small-scale model with nine different excitations, and the matrix exponential of the small-scale model just need to be computed once. For the remaining small-scale models, because the boundary conditions of each smallscale model may be different, the responses of these models need to be 221

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

conduction, computation of the solutions of the basic finite element models can be transformed into computation of the solutions of the small-scale models b(kx , ky) (k x = 1,2, ⋯, Nx ; k y = 1,2, ⋯, Ny ) . Because we assume that excitations can only affect m unit cells in any direction within one time step, each small-scale model b(kx , ky) is composed of unit cells. Furthermore, when (2m + 1) × (2m + 1) m + 1 ≤ k x ≤ Nx − m and m + 1 ≤ k y ≤ Ny − m , arbitrary homogenous boundary conditions can be applied on the boundaries of the smallscale model b(kx , ky) , while for the other small-scale models, parts of the boundary conditions are the same as those of the basic finite element model, and other boundary conditions can be set as arbitrary homogenous boundary conditions. Finally, since arbitrary homogenous boundary conditions can be applied on the boundaries of each small-scale model b(kx , ky) (m + 1 ≤ k x ≤ Nx − m and m + 1 ≤ k y ≤ Ny − m ), we can apply the same homogenous boundary conditions for all these models. Then, according to the periodic property of the structure, the (Nx − 2m)(Ny − 2m) small-scale models are the same except that their excitations may be different. Therefore, computation of the responses of the (Nx − 2m)(Ny − 2m) small-scale models can be computation of computing the responses of one small-scale model with (Nx − 2m)(Ny − 2m) different excitations, and the matrix exponential of the small-scale model just needs to be computed once. For the smallscale model b(kx , ky) (1 ≤ k x ≤ m and 1 ≤ k y ≤ m ), if the same homogenous boundary conditions are applied on the top and right boundaries of the small-scale model b(kx , ky) , then the m2 small-scale models are the same except that their excitations may be different. Therefore,

calculated individually. In this subsection, according to the superposition principle and the physical features of the transient heat conduction and the periodic property of the structure, for a 2D periodic structure with 5 × 5 unit cells, we have shown that computation of the temperature response of the original finite element models can be transformed into computation of the temperature responses of the small-scale models. Obviously, the idea used for this example can be generalized directly to the general 2D periodic structures without any difficulty, which is summarized in the following subsection.

5.2. The method for general 2D periodic structures In this subsection, the same idea proposed in subsection 5.1 is generalized to the general 2D periodic structures, which leads to an efficient and accurate method for solving the transient heat conduction in 2D periodic structures. First, based on the superposition principle of linear systems, the solution of the original finite element model can be obtained by the superposition of the solutions of the Nx × Ny basic finite element models B (kx , ky) (k x = 1,2, ⋯, Nx ; k y = 1,2, ⋯, Ny ) . The solution domain of the basic finite element model B (kx , ky) is the same as that of the original finite element model, but the excitations of the basic finite element (k , k ) model B (kx , ky) are the basic external excitation r α x y and the basic in(kx , ky ) . itial temperature T k Second, according to the physical feature of transient heat

Fig. 7. The small-scale model consisting of (2m + 1) × (2m + 1) unit cells and its interfaces #1 and #m + 1.

222

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

be used to the final time step. The detailed strategy is as follows. First, the desired parameter m and an initial time step η1 are given, and the heat capacity matrix C and the heat conduction matrix K of a small-scale model are formed; then, the matrix H of the small-scale model is computed. Second, the parameters N and q can be determined based on the matrix H and the initial time step η1 [87–89]. Then, according to Eq. (16), the matrix exponential S0 = I + R 0 corresponding to the time step η1/2N is computed. Next, we determine whether the responses caused by the excitations can spread over m unit cells. If they can, then the initial time step η1 is too large, and thus the programme should be stopped and restarted by adjusting the value of η1; otherwise, perform the loop (17). In each cycle of the loop (17), we determine whether the responses caused by the excitations can spread over m unit cells according to the matrix R i . When they can, exit the loop after storing the matrix R i − 1 and the time step η1⋅2(i − N − 1) ; otherwise, perform the next cycle of the loop (17). In the above procedure, the key step is to decide whether the responses caused by the excitations within a unit cell can spread over m unit cells, which can be accomplished as follows. First, a small parameter ε is chosen, for example, ε = 10−16 , and s#1 and s# m + 1 respectively denote the nodes on interface #1 and #m + 1 (see Fig. 7). For the i−th cycle of the loop (17), the time step is η1⋅2(i − N ) , and the corresponding matrix exponential is Si = I + R i . According to the physical meaning of the matrix exponential, we can use Si (s# m + 1, s#1) to determine whether the responses caused by the excitations can spread over m unit cells, i.e., Si (s# m + 1, s#1) ∞ ≥ ε means that the response has spread over m unit cells, and otherwise it has not yet spread over m unit cells. According to the above strategy, once a desired parameter m is given, a reasonable time step η and the matrix exponential corresponding to the small-scale model can be obtained. The choice of the parameter m is discussed as follows. On the one hand, according to the discussion in Section 4, by using a smaller scale model, we can make full use of the periodic property of the structure and the sparsity of the matrix exponential. On the other hand, if the size of the small-scale model is larger, the automatically selected time step becomes larger. Thus, when the time interval of the simulation is fixed, the total number of steps becomes less, while the CPU time for each time step becomes larger. Based on the above two points, to select a reasonable time step and make full use of the periodic property of the structure and the sparsity of the matrix exponential, the parameter m should be set as a small value.

computation of the responses of the m2 small-scale models can be transformed into computation of the responses of one small-scale model with m2 different excitations, and the matrix exponential of the smallscale model just needs to be computed once. For other small-scale models, similar analysis can be applied. Eventually, only 2Nx + 2Ny − 8m + 5 matrix exponentials corresponding to the smallscale models need to be computed to solve the transient heat conduction problem of a 2D periodic structure with Nx × Ny unit cells. Remarkably, the boundary conditions applied on the 2D periodic structure are usually periodic in actual engineering applications, i.e., the boundary conditions that are applied on each unit cell on an edge of the 2D periodic structures are all same. As shown in subsection 5.1, arbitrary homogenous boundary conditions can be applied on some boundaries of the small-scale models b(kx , ky) . Using this property and the fact that the boundary conditions are periodic, the temperature response of the original periodic structure can be obtained by computing the temperature responses of one small-scale model with different excitations. Therefore, the matrix exponential of the small-scale model for this case only needs to be computed once. In summary, for a 2D periodic structure with general boundary conditions or periodic boundary conditions, only 2Nx + 2Ny − 8m + 5 or one matrix exponential corresponding to the small-scale model with a few unit cells needs to be computed to solve the transient heat conduction of the original 2D periodic structure. When the scale of the original 2D periodic structure is very large and the value of m is small, the size of the matrix exponential of the small-scale model is much less than that of the original periodic structure. Therefore, the proposed algorithm can significantly improve the computational efficiency and reduce the memory cost. 6. Determination of the parameter m and the time step η Based on the analysis of Section 5, for the proposed method, the temperature responses can be computed by using matrix exponentials of some small-scale models with (m + 1) × (m + 1) unit cells instead of that of the original periodic structure. Once the parameter m and the time step η are determined, the matrix exponential of each small-scale model can be computed easily according to PIM introduced in Section 3. However, we suppose that the excitations can only affect m unit cells in any direction within one time step η in the previous analysis, and the explicit relationship between the parameter m and the time step η cannot be known in advance. Thus, the following strategy is proposed to determine the parameter m and the time step η. A desired parameter m and an initial time step η1 are given such that it is believed that the excitations of a unit cell can affect m unit cells in any directions within one time step η1. If the shortest time for the excitations to affect m unit cells is found to be η2 , then η = min(η1, η2) can

7. Numerical examples Example 1. A 2D periodic structure composed of N × N unit cells is considered in this example. The unit cell is a square with a square hole Fig. 8. (a) The unit cell for Example 1; (b) the FEM mesh of the unit cell.

223

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

proposed method. When t = 3200 s , the temperature contours computed by the Crank-Nicholson method with time steps of 16, 8, and 0.01 s and the proposed method are shown in Fig. 9(a)-(d), respectively. Fig. 9 shows that the temperatures computed by the proposed method with a time step of 16 s and by the Crank-Nicholson method with a time step of 0.01 s are almost identical, which verifies the correctness of the proposed method. Alternately, Fig. 9 also shows that the results computed by the Crank-Nicholson method with the same time step as the proposed method are incorrect and that, even when the time step of the Crank-Nicholson method is twice as small as that of the proposed method, the differences between the results obtained from the CrankNicholson method and the reference solutions cannot be ignored. Based on the above analysis, it can be concluded that the proposed method can give very accurate solutions using a larger time step, while only with a relatively small time step can the Crank-Nicholson method give more reasonable solutions. To further illustrate the accuracy of the proposed method, the relative errors of the results obtained from the proposed method and the Crank-Nicholson method are defined as

and consists of two different materials denoted as #1 and #2. The sizes of the unit cell are shown in Fig. 8(a). The heat conduction coefficients for the two materials are k1 = 25 Wm−1K−1 and k2 = 0.2 Wm−1K−1, and the volumetric heat capacities are ρ1 c1 = 3 × 106 Jm−3K−1 and ρ2 c2 = 1 × 106 Jm−3K−1, respectively. The unit cell is discretized using the FEM with three-node triangle linear elements, and the mesh is shown in Fig. 8(b). The initial condition of the 2D periodic structure is

T (x , y, 0) = 30 °C

(28)

The heat exchanges are modelled by a Robin condition at each boundary surface of the periodic structure with the external temperature Text = 200 °C and a heat transfer coefficient h = 500 Wm−2K−1, i.e., ∂T

∂T

− k ∂x (0, y, t ) = −k ∂y (x , 0, t ) = h (Text − T ) Wm−2 ∂T

∂T

k ∂x (0.1N , y, t ) = k ∂y (x , 0.1N , t ) = h (Text − T ) Wm−2

(29)

The boundary condition for each square hole is kept insulated. The volumetric heat source Q (x , y, t ) is imposed at the centroid of the periodic structure, i.e.,

eR = max[ (T − TREF )/ TREF ]

Q (x , y, t ) = δ (x − 0.05N , y − 0.05N ) f (t ) 10 4 (2

f (t ) =

⎧ ⎪ ⎨ ⎪0 ⎩

− sin(πt /16))

Wm−2,

where TREF represents the reference solutions of all nodes at a given time using the Crank-Nicholson method with a time step of 0.01 s, and T represents the temperature responses of all nodes at the same time computed using the Crank-Nicholson method with time steps of 16, 8 and 4 s and the proposed method with a time step of 16 s. When N = 20 , the relative errors at different times are given in Table 1, and the values corresponding to the Crank-Nicholson method with time steps of 16, 8 and 4 s are more than 49%, 17% and 4%, while the values of the proposed method with a time step of 16 s are approximately 10−7 . This shows that the precision of the Crank-Nicholson method with a larger time step is very low, while the proposed method can give very accurate solutions by using a larger time step. To illustrate the efficiency of the proposed method, we define a

16(n − 1) ≤ t < 16n (n = 1,3,5, ⋯) otherwise

(31)

(30)

The transient temperature responses are computed using the proposed method and the Crank-Nicholson method. The interval of integration is from 0 to 9600 s. The small-scale model is composed of 3 × 3 unit cells for the proposed method, and the time step selected automatically by the proposed method is 16 s. Three different time steps of 16, 8 and 4 s are used for the Crank-Nicholson method, and the numerical result computed using the Crank-Nicholson method with a time step of 0.01 s is regarded as the reference solution. The case of N = 10 is used to prove the high accuracy of the

Fig. 9. Comparisons of the temperature contours obtained by the Crank-Nicholson method with different time steps and the proposed method when N = 10 and t = 3200 s for Example 1: (a)–(c) respectively represent the temperature contours computed by the Crank-Nicholson method with time steps of 16, 8 and 0.01 s, and (d) represents the temperature contours computed by the proposed method with a time step of 16 s.

224

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

Table 1 Relative errors (%) at different times for the proposed and Crank-Nicholson methods for Examples 1 and 2. Time

32 s 4800 s 9600 s

Example 1

Time

Proposed method (η )

Crank-Nicholson method (η )

16s

16 s

8s

4s

8.8E-6 1.5E-5 1.5E-5

146.3 50.0 49.5

31.5 17.8 17.0

8.9 4.7 4.5

2s 1600 s 3200 s

Example 2 Proposed method (η )

Crank-Nicholson method (η )

2s

2s

1s

0.5 s

6.7E-4 1.5E-3 1.5E-3

33.5 74.6 73.2

8.6 16.2 15.9

2.2 3.9 3.8

relative time tR as

tR = tCN / tFPIM

(32)

where tCN and tFPIM denote the CPU times for the Crank-Nicholson method and the proposed method. The interval of integration is [0,9600s]. The direct method and the preconditioning conjugate gradient (PCG) method are used to solve the system of linear equations obtained by the Crank-Nicholson method. For the direct method, using Cholesky decomposition, the symmetric positive-definite coefficient matrix of the system of linear equations is decomposed into a lower triangular matrix, and then the system of linear equations is solved using a backsubstitution approach. For the PCG method, the incomplete Cholesky decomposition is applied to the symmetric positive-definite matrix of the system of linear equations, and then the system of linear equations is solved using the conjugate gradient method. In this example, both the absolute and relative tolerances for the PCG method are set as 10−6 . The curves of relative time tR corresponding to the direct and PCG methods with different time steps are shown in Fig. 10, in which the

Fig. 10. The relative times for Example 1 corresponding to the direct method with time steps of 16 s (diamond), 8 s (hexagon), and 4 s (circle) and the PCG method with time steps of 16 s (pentagram), 8 s (triangle), and 4 s (square).

Fig. 11. The four types of unit cells of lattice periodic structures.

225

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

Fig. 12. The FEM meshes of the four types of unit cells.

δ = 0.01 m α = 2δ

lines with diamonds, hexagons and circles are the relative time tR corresponding to the direct method with time steps of 16, 8 and 4 s, respectively, and the lines with pentagrams, triangles and squares are the relative time tR corresponding to the PCG method with time steps of 16, 8 and 4 s, respectively. When the direct method is used to solve the system of linear equations obtained by the Crank-Nicholson method, Fig. 10 shows that the relative time tR increases along with the increase of N , and for N = 90 , the efficiency of the proposed method is approximately 24, 44 and 74 times better than that of the Crank-Nicholson method with time steps of 16, 8 and 4 s; for N > 90 , the direct method cannot be performed due to computer memory limitations. In addition, when the PCG method is used to solve the system of linear equations obtained by the Crank-Nicholson method, Fig. 10 shows that, for a larger value of N , the efficiency of the proposed method is approximately 11, 17 and 28 times better than that of the Crank-Nicholson method with time steps of 16, 8 and 4 s. Note that, when N = 200 , the periodic structure is composed of 200 × 200 unit cells and has 5884001 nodes and 10880000 elements. Based on Table 1 and Fig. 10, it can be observed that the proposed method is very accurate and efficient, while the Crank-Nicholson method with the same or smaller time steps cannot achieve such precision and consumes significant of CPU resources. Therefore, in terms of precision and efficiency, the proposed method has a huge advantage. Example 2. In this example, the 2D transient heat conduction problems of lattice periodic structures are considered. The four types of unit cells of the periodic structures are shown in Fig. 11(a)-(d). The unit cells can be regarded as a square plate that is removed from some solid material with the same volume ratio and different shapes. The length of each unit cell is 0.1 m, and other parameters of the unit cells shown in Fig. 11 are

2

β=

(4.2 + 1.2 3 ) − (1824 3 + 1152)(0.1 − (1 + 2 ) δ ) + (3.84 3 − 7.68)

γ = (0.5 +

114 + 24 3

2 /2) δ

(33)

The heat conduction coefficient and the volumetric heat capacity are k = 200 Wm−1K−1 and ρc = 2.5 × 106 Jm−3K−1, respectively. The four types of unit cells are discretized using the FEM with three-node triangle linear elements (see Fig. 12(a)-(d)). The 2D periodic structure is composed of N × N unit cells, and the initial condition of the periodic structure is

T (x , y, 0) = 30 °C

(34)

The temperature is prescribed at the boundary surface y = 0.1N , i.e.,

T (x , 0.1N , t ) = 30 °C

(35)

An adiabatic condition is assumed at the boundary surface y = 0 , i.e.,

k

∂T (x , 0, t ) = 0 ∂y

(36)

The heat exchanges are modelled by a Robin condition at the boundary surface x = 0.1N with the external temperature Text = 30 °C and a heat transfer coefficient h = 500 Wm−2K−1, i.e.,

k

∂T (0.1N , y, t ) = h (Text − T ) Wm−2 ∂x

The boundary condition at the boundary surface x = 0 is 226

(37)

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

Fig. 13. The temperature contours of the four types of periodic structures at t = 3200 s .

Q (x , y, t ) =

⎧ q (t ), when (N − 1)/20 ≤ x ≤ (N + 1)/20 and (N − 1) ⎪ /20 ≤ y ≤ (N + 1)/20 ⎨ ⎪ 0, otherwise ⎩

107 Wm−2, 2(n − 1) ≤ t < 2n q (t ) = ⎧ ⎨ otherwise ⎩0

(n = 1,3,5, ⋯) (39)

The transient heat conduction problems of the four types of periodic structures are solved using the proposed method. The small-scale model is composed of 3 × 3 unit cells, and the time step selected automatically by the proposed method is 2 s. The interval of integration is from 0 to 3200 s. When the parameter N is set as 11, i.e., there are 121 unit cells, the temperature contours corresponding to the four types of periodic structures at t = 3200 s are as shown in Fig. 13(a)-(d). It can be seen that there are significant differences between the temperature contours for the four different periodic structures, which means that, when the volume ratio is the same, the thermal conductivity strongly depends on the geometric shape of the lattice [28,90]. To further illustrate the accuracy and efficiency of the proposed method, we focus on the periodic structure consisting of the unit cell shown in Fig. 11(a). The temperature responses of the periodic structures are calculated with the proposed method and the Crank-Nicholson method. The size of the small-scale model, the time step of the proposed method, and the integration interval are all the same as mentioned above. Three time steps of 2, 1 and 0.5 s are used for the Crank-Nicholson method. The numerical results computed by the Crank-Nicholson method with a time step of 0.01 s are used as the reference solution. When N = 11, the relative errors of the results obtained from the proposed method and the Crank-Nicholson method with the different time steps of 2, 1 and 0.5 s are given in Table 1. We can see that, as time goes on, the relative errors of the results obtained from the Crank-

Fig. 14. The relative times for Example 2 corresponding to the direct method with time steps of 2 s (diamond), 1 s (hexagon), and 0.5 s (circle) and the PCG method with time steps of 2 s (pentagram), 1 s (triangle), and 0.5 s (square).

−

∂T k ∂x

⎧ f (t ) + h (Text − T ) Wm−2, when (N − 3)/20 ≤ y ⎪ (0, y, t ) = ≤ (N + 3)/20 ⎨ ⎪ h (Text − T ) Wm−2, otherwise ⎩ f (t ) = 105 [2 − 1.5 × sin(πt /2) ] Wm−2 (38)

The boundary condition for each internal pore is kept insulated. The following volumetric heat source Q (x , y, t ) is imposed on the periodic structure:

227

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

Fig. 15. (a) The unit cell and its FEM mesh for Example 3; (b) in-plane size of the unit cell.

Fig. 16. The temperature contours of the front surface when N = 21 and t = 16000 s for Example 3: (a)–(c) respectively represent the temperature contours computed by the Crank-Nicholson method with time steps of 8, 4 and 2 s, and (d) represents the temperature contours computed by the proposed method with a time step of 8 s.

lines with diamonds, hexagons and circles correspond to the direct method with time steps of 2, 1 and 0.5 s, and the solid lines with pentagrams, triangles and squares correspond to the PCG method with time steps of 2, 1 and 0.5 s, respectively. It can be seen that when the direct method is used to solve the system of linear equations obtained by the Crank-Nicholson method, the relative time tR increases with the increase of N . When N = 41, the efficiency of the proposed method is approximately 12, 22 and 63 times better than that of the Crank-Nicholson method with time steps of 2, 1 and 0.5 s, respectively. When N > 41, the direct method cannot be performed due to the computer memory limitations. In addition, when the PCG method is used to solve the system of linear equations obtained by the Crank-Nicholson method, for a larger N , the efficiency of the proposed method is approximately 13, 20 and 32 times better than that of the Crank-

Nicholson method with time steps of 2, 1 and 0.5 s are approximately 74%, 16% and 4%, while the relative errors of the results obtained from the proposed method with a time step of 2 s are less than 10−4 . Therefore, it can be concluded that the proposed method, even with a larger time step, is very accurate, while the Crank-Nicholson method with the same step cannot give satisfactory results. Next, Eq. (32) is used to compare the efficiency of the proposed method and the Crank-Nicholson method with different time steps. The interval of integration is [0,3200 s]. The direct method and the Preconditioning Conjugate Gradient (PCG) method are used to solve the system of linear equations obtained by the Crank-Nicholson method, and both the absolute and relative tolerances for the PCG method are set as 10−6 . Fig. 14 shows the curves of relative time tR corresponding to the direct and PCG methods with different time steps, in which the solid 228

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

solve this problem. The integration interval is [0, 16000 s]. The smallscale model is composed of 3 × 3 unit cells for the proposed method, and the time step selected automatically by the proposed method is 8 s. Three different time steps of 8, 4 and 2 s are used for the CrankNicholson method. The PCG method is used to solve the system of linear equations obtained by the Crank-Nicholson method, and both the absolute and relative tolerances for the PCG method are set as 10−6 . When t = 16000 s , the temperature contours of the front surface, i.e., z = 0.05 m , computed by the Crank-Nicholson method with time steps of 8, 4, and 2 s and the proposed method are shown in Fig. 16(a)-(d), respectively. Fig. 16 shows that if the Crank-Nicholson method uses the same time step as the proposed method, i.e., 8 s, the results are incorrect; when the time step of the Crank-Nicholson method is two times smaller than that of the proposed method, the differences between the temperature contours obtained from the Crank-Nicholson method and the proposed method are still obvious; when the time step of the CrankNicholson method is 4 times as small as the time step of the proposed method, the results given by the Crank-Nicholson method are almost the same as those given by the proposed method. Based on the above analysis, it can be concluded that only with the relatively smaller time step can the Crank-Nicholson method give more reasonable solutions, while the proposed method can give very accurate solutions using a larger time step. To further illustrate the accuracy of the proposed method, an error e is defined as

Fig. 17. Errors of the results obtained from the Crank-Nicholson method with time steps of 8 s (diamond), 4 s (square) and 2 s (triangle) for Example 3.

Table 2 Comparison of CPU times for example 3. Methods

Proposed method

Crank-Nicholson method

Time step (s) Steps CPU time (s)

8 2000 345

8 2000 1847

4 4000 2894

2 8000 4733

e = max[ (T − TFPIM )/ TFPIM ]

Nicholson method with the three time steps. Note that, when N = 201, the periodic structure is composed of 201 × 201 unit cells and has 18430897 nodes and 30704760 elements. Table 1 and Fig. 14 again show that, to obtain reasonable results, a smaller time step should be used for the Crank-Nicholson method, while a much larger time step can be used for the proposed method. Therefore, to achieve the same precision, the proposed method is much more efficient than the Crank-Nicholson method. Example 3. The 3D transient heat conduction of a periodic structure is considered. The unit cell of the periodic structure is shown in Fig. 15(a). The thickness of the unit cell is 0.1 m in the z direction, and the plane sizes are shown in Fig. 15(b). The heat conductivity coefficient and the volumetric heat capacity of the material are k = 80 Wm−1K−1 and ρc = 3.537 × 106 Jm−3K−1, respectively. The unit cell is discretized using the FEM, and the mesh is shown in Fig. 15(a). The periodic structure is composed of 21 × 21 unit cells and the computational domain is from 0 to 2.1 m in direction x, from 0 to 2.1 m in direction y, and from −0.05 to 0.05 m in direction z. The discretized periodic structure has 550589 nodes and 2050650 elements. The initial condition of the periodic structure is

T (x , y, z , 0) = 30 °C

where TFPIM denotes the temperature responses of all nodes at a given time computed using the proposed method, and T denotes the temperature responses of all nodes at the same time computed using the Crank-Nicholson method with the alternative time steps. The errors corresponding to the results obtained from the Crank-Nicholson method are shown in Fig. 17, in which the solid lines with diamonds, squares and triangles correspond to the time steps of 8, 4 and 2 s, respectively. Fig. 17 shows that, if the time step of the Crank-Nicholson method is smaller, the results obtained from the Crank-Nicholson method are closer to those obtained from the proposed method. Hence, the correctness of the proposed method is proved again. To illustrate the efficiency of the proposed method, Table 2 gives the CPU times used by the proposed and Crank-Nicholson methods, which shows that the efficiency of the proposed method is approximately 5, 8 and 14 times better than that of the Crank-Nicholson method with the three time steps. This example again shows that the proposed method is more efficient and accurate than the Crank-Nicholson method. 8. Conclusions

(40)

An efficient and accurate method is developed to analyse the transient heat conduction in the 2D periodic structure. Based on the superposition principle in linear systems, the physical feature of transient heat conduction and the periodic property of the structure, computation of the temperature responses of the original large-scale 2D periodic structure is transformed into computation of the temperature responses of the small-scale models, which can significantly improve computational efficiency and reduce the memory requirement. Because the above transformation is accurate and the responses of the small-scale models are obtained using PIM, the proposed method inherits the accuracy and stability of the PIM, which can yield a highly accurate solution with a larger time step. The following conclusions can be drawn based on the three numerical examples. With the larger time steps, the efficiency of the proposed method is approximately more than 10 and 5 times better than the Crank-Nicholson method with the direct solver and the PCG solver, respectively. In terms of precision, even if the Crank-Nicholson method uses a time step which is four times as small as that of the proposed method, its precision is less than one-thousandth of

The temperature is prescribed at the boundary surfaces x = 2.1 m and y = 2.1 m , i.e.,

T (x , 2.1, z , t ) = T (2.1, y, z , t ) = 100 °C

(43)

(41)

On the other boundary surfaces, the following heat flux q is imposed on the structure

⎧ f (t ), when x = 0, 1 ≤ y ≤ 1.1 and − 0.05 ≤ z ≤ 0.05 q (x , y, z , t ) = 5f (t ), when 1 ≤ x ≤ 1.1 and 1 ≤ y ≤ 1.1 and z = 0.05 ⎨ otherwise ⎩ 0, ⎧105 [2 − sin(πt /16) ] Wm−2, 16(n − 1) ≤ t < 16n (n ⎪ ⎪ = 1,3,5, ⋯) f (t ) = ⎨ 16(n − 1) ≤ t < 16n (n otherwise ⎪ ⎪ = 1,3,5, ⋯) ⎩ (42) The proposed method and the Crank-Nicholson method are used to 229

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

that of the proposed method. Thus, to compare with the CrankNicholson method, the proposed method can solve the transient heat conduction problems of 2D periodic structures more accurately and effectively.

pores using a homogenization method, Heat Mass Tran 52 (2016) 2213–2217. [29] Y. Asakuma, T. Yamamoto, Thermal analysis of resin composites with ellipsoidal filler considering thermal boundary resistance, J Therm Sci 25 (2016) 424–430. [30] Y. Asakuma, Y. Kanazawa, T. Yamamoto, Thermal radiation analysis of packed bed by a homogenization method, Int J Heat Mass Transfer 73 (2014) 97–102. [31] F.M. Zhai, L.Q. Cao, Multiscale asymptotic method for heat transfer equations in lattice-type structures, Int J Numer Anal Mod 6 (2009) 232–255. [32] K. Terada, M. Kurumatani, T. Ushida, N. Kikuchi, A method of two-scale thermomechanical analysis for porous solids with micro-scale heat transfer, Comput Mech 46 (2010) 269–285. [33] Q. Ma, J.Z. Cui, Second-order two-scale analysis method for the heat conductive problem with radiation boundary condition in periodical porous domain, Commun Comput Phys 14 (2013) 1027–1057. [34] Z.Q. Yang, J.Z. Cui, Z.Q. Wang, Y. Zhang, Multiscale computational method for nonstationary integrated heat transfer problem in periodic porous materials, Numer. Meth. PDEs 32 (2016) 510–530. [35] Z.Q. Yang, J.Z. Cui, Y. Sun, J.G. Ge, Multiscale computation for transient heat conduction problem with radiation boundary condition in porous materials, Finite Elem Anal Des 102 (2015) 7–18. [36] Z.Q. Yang, J.Z. Cui, Y. Sun, Transient heat conduction problem with radiation boundary condition of statistically inhomogeneous materials by second-order twoscale method, Int J Heat Mass Transfer 100 (2016) 362–377. [37] Z.Q. Yang, J.Z. Cui, Y.F. Nie, Q. Ma, The second-order two-scale method for heat transfer performances of periodic porous materials with interior surface radiation, Comput Model Eng Sci: Comput Model Eng Sci 88 (2012) 419–442. [38] Z.Q. Yang, J.Z. Cui, S. Zhou, Thermo-mechanical analysis of periodic porous materials with microscale heat transfer by multiscale asymptotic expansion method, Int J Heat Mass Transfer 92 (2015) 904–919. [39] A. Dasgupta, R.K. Agarwal, Orthotropic thermal conductivity of plain-weave fabric composites using a homogenization technique, J Compos Mater 26 (1992) 2736–2758. [40] A. Alzina, E. Toussaint, A. Béakou, B. Skoczen, Multiscale modelling of thermal conductivity in composite materials for cryogenic structures, Compos Struct 74 (2006) 175–185. [41] I. Ozdemir, W. Brekelmans, M. Geers, Computational homogenization for the thermo-mechanical analysis of heterogeneous solids, Comput Meth Appl Mech Eng 112 (2008) 602–613. [42] V.G. Kouznetsova, M.G.D. Geers, W.A.M. Brekelmans, Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy, Comput Meth Appl Mech Eng 193 (2004) 5525–5550. [43] B.R. Wang, J.T. Liu, S.T. Gu, Q.C. He, Numerical evaluation of the effective conductivities of composites with interfacial weak and strong discontinuities, Int J Therm Sci 93 (2015) 1–20. [44] X.Y. Cui, Z.C. Li, H. Feng, S.Z. Feng, Steady and transient heat transfer analysis using a stable node-based smoothed finite element method, Int J Therm Sci 110 (2016) 12–25. [45] A. Bensoussan, J.L. Lions, G. Papanicolaou, Asymptotic analysis for periodic structures, North-Holland, Amsterdam, 1978. [46] J.L. Lions, Some methods in the mathematical analysis of systems and their control, Science Press, Beijing, 1981. [47] P. Ylä-Oijala, S. Järvenpää, Iterative solution of high-order boundary element method for acoustic impedance boundary value problems, J Sound Vib 291 (2006) 824–843. [48] S. Järvenpää, M. Taskinen, P. Ylä-Oijala, Singularity extraction technique for integral equation methods with higher order basis functions on plane triangles and tetrahedra, Int J Numer Meth Eng 58 (2003) 1149–1165. [49] R.W. Freund, Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices, SIAM J Sci Statist Comput 13 (1992) 425–448. [50] A. Alia, H. Sadok, M. Souli, CMRH method as iterative solver for boundary element acoustic systems, Eng Anal Bound Elem 36 (2012) 346–350. [51] H. Xiao, Z. Chen, Numerical experiments of preconditioned Krylov subspace methods solving the dense non-symmetric systems arising from BEM, Eng Anal Bound Elem 31 (2007) 1013–1023. [52] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J Sci Stat Comput 7 (1986) 856–869. [53] C. Serra, A. Tadeu, N. Simões, Iterative simulation of 3D heat diffusion in a medium with multiple cracks, Eng Anal Bound Elem 41 (2014) 10–17. [54] J.M. Winget, T.J.R. Hughes, Solution algorithms for nonlinear transient heat conduction analysis employing element-by-element iterative strategies, Comput Meth Appl Mech Eng 52 (1985) 711–815. [55] P. Sun, Z. Luo, Y. Zhou, Some reduced finite difference schemes based on a proper orthogonal decomposition technique for parabolic equations, Appl Numer Math 60 (2010) 154–164. [56] R.A. Biaecki, A.J. Kassab, A. Fic, Proper orthogonal decomposition and modal analysis for acceleration of transient FEM thermal analysis, Int J Numer Meth Eng 62 (2005) 774–797. [57] A. Fic, R.A. Białecki, A.J. Kassab, Solving transient nonlinear heat conduction problems by proper orthogonal decomposition and the finite-element method, Numer Heat Transfer B 48 (2005) 103–124. [58] P. Ding, X.H. Wu, Y.L. He, W.Q. Tao, A fast and efficient method for predicting fluid flow and heat transfer problems, J Heat Tran 130 (2008) 171–181. [59] X. Zhang, X. Hui, A fast meshless method based on proper orthogonal decomposition for the transient heat conduction problems, Int J Heat Mass Transfer 84 (2015) 729–739. [60] X.H. Zhang, J. Ouyang, L. Zhang, Matrix free meshless method for transient heat conduction problems, Int J Heat Mass Transfer 52 (2009) 2161–2165.

Acknowledgements The authors would like to thank the editor and three anonymous reviewers for their valuable comments and suggestions improved the overall quality of the paper. The authors are grateful for the financial support of the Natural Science Foundation of China (No. 11572076) and the 973 program (No. 2014CB049000). References [1] L. Desgrosseilliers, D. Groulx, M.A. White, Heat conduction in laminate multilayer bodies with applied finite heat source, Int J Therm Sci 72 (2013) 47–59. [2] P.B. Deshpande, J.R. Couper, Thermal conductivity of two-phase systems, J Heat Tran 94 (1972). [3] H.S. Carslaw, J.C. Jaeger, H. Feshbach, Conduction of heat in solids, Oxford University Press, Oxford, 1959. [4] G. Horvay, R. Mani, M.A. Veluswami, G.E. Zinsmeister, Transient heat conduction in laminated composites, J Heat Tran 95 (1973) 309–316. [5] A.M. Manaker, G. Horvay, Thermal response in laminated composites, Z Angew Math Mech 55 (1975) 503–513. [6] A.H. Nayfeh, A continuum mixture theory of heat conduction in laminated composites, J Appl Mech 42 (1975) 399–404. [7] J.L. Auriault, Effective macroscopic description for heat conduction in periodic composites, Int J Heat Mass Transfer 26 (1983) 861–869. [8] S.J. Matysiak, C. Woźniak, On the modelling of heat conduction problem in laminated bodies, Acta Mech 65 (1986) 223–238. [9] S.J. Matysiak, On certain problems of heat conduction in periodic composites, Z Angew Math Mech 71 (1991) 524–528. [10] S.J. Matysiak, O.M. Ukhanska, On two-dimensional heat conduction problem of a periodic stratified convective half-space with a moving heat source, Int Commun Heat Mass Transfer 24 (1997) 129–138. [11] A.K. Noor, L.H. Tenek, Steady-state nonlinear heat transfer in multilayered composite panels, J Eng Mech 118 (1992) 539–542. [12] C. Woźaniak, Z.F. Baczyński, M. Woźniak, Modelling of nonstationary heat conduction problems in micro-periodic composites, Z Angew Math Mech 76 (1996) 285–287. [13] S.J. Matysiak, A.A. Yevtushenko, Mixed nonstationary problem of heat conduction for a microperiodic two-layered half-space, Int Commun Heat Mass Transfer 34 (2007) 1101–1107. [14] J. Ignaczak, Z.F. Baczyński, On a refined heat conduction theory for microperiodic layered solids, J Therm Stresses 20 (1997) 749–771. [15] H. Bufler, Stationary heat conduction in a macro- and microperiodically layered solid, Arch Appl Mech 70 (2000) 103–114. [16] G. Laschet, Homogenization of the thermal properties of transpiration cooled multilayer plates, Comput Meth Appl Mech Eng 191 (2002) 4535–4554. [17] R. Kulchytsky-Zhyhailo, S.J. Matysiak, On heat conduction problem in a semi-infinite periodically laminated layer, Int Commun Heat Mass Transfer 32 (2005) 123–132. [18] R. Kulchytsky-Zhyhailo, S.J. Matysiak, On some heat conduction problem in a periodically two-layered body. Comparative results, Int Commun Heat Mass Transfer 32 (2005) 332–340. [19] R. Kulchytsky-Zhyhailo, S.J. Matysiak, On temperature distributions in a semi-infinite periodically stratified layer, Bull. Pol. Ac.: Tech 54 (2006) 45–49. [20] G.F. Abdelal, M.M. Abdalla, Z. Gürdal, Transient heat conduction of variable stiffness composite laminate, J Therm Stresses 33 (2010) 121–134. [21] A.H. Muliana, J.S. Kim, A two-scale homogenization framework for nonlinear effective thermal conductivity of laminated composites, Acta Mech 212 (2010) 319–347. [22] P.W. Chung, K.K. Tamma, R.R. Namburu, Homogenization of temperature-dependent thermal conductivity in composite materials, J Thermophys Heat Transfer 15 (2012) 10–17. [23] A. Matine, N. Boyard, P. Cartraud, G. Legrain, Y. Jarny, Modeling of thermophysical properties in heterogeneous periodic media according to a multi-scale approach: effective conductivity tensor and edge effects, Int J Heat Mass Transfer 62 (2013) 586–603. [24] A. Matine, N. Boyard, G. Legrain, Y. Jarny, P. Cartraud, Transient heat conduction within periodic heterogeneous media: a space-time homogenization approach, Int J Therm Sci 92 (2015) 217–229. [25] X. Wang, X.l. Duan, Y. Gao, Multiscale asymptotic analysis and parallel algorithm of parabolic equation in composite materials, Math Proble Eng 2014 (2014). [26] E. Monteiro, J. Yvonnet, Q.C. He, Computational homogenization for nonlinear conduction in heterogeneous materials using model reduction, Comp Mater Sci 42 (2008) 704–712. [27] S. Liu, Y. Zhang, Multi-scale analysis method for thermal conductivity of porous material with radiation, Multidiscip Model Mater Struct 2 (2006) 327–344. [28] Y. Asakuma, T. Yamamoto, Thermal analysis of porous medium with ellipsoidal

230

International Journal of Thermal Sciences 127 (2018) 213–231

Q. Gao, H.C. Cui

[78] J.H. Lin, W.X. Zhong, W.S. Zhang, D.K. Sun, High efficiency computation of the variances of structural evolutionary random responses, Shock Vib 7 (2015) 209–216. [79] S.J. Tan, H.J. Peng, W.Y. Zhou, Z.G. Wu, A novel extended precise integration method based on Fourier series expansion for the H2-norm of linear time-varying periodic systems, Int J Control (2016) 1–25. [80] Z. Cai, Y. Zhang, X. Li, B. Jin, Research on the vibration control algorithms for large space truss structure based on precise integration method, Int J Comput Meth 13 (2016) 1641021. [81] T.C. Fung, A precise time-step integration method by step-response and impulsiveresponse matrices for dynamic problems, Int J Numer Meth Eng 40 (1997) 4501–4527. [82] J.H. Lin, W.P. Shen, F.W. Williams, A high precision direct integration scheme for structures subjected to transient dynamic loading, Comput Mech 56 (1998) 113–120. [83] W.X. Zhong, F.W. Williams, A precise time step integration method. Proceedings of the intistution of mechanical engineers part C, J Mech Eng Sci 208 (1994) 427–430. [84] W.X. Zhong, J.N. Zhu, X.X. Zhong, On a new time integration method for solving time dependent partial differential equations, Comput Meth Appl Mech Eng 130 (1996) 163–178. [85] A. Tadeu, N. Simões, Three-dimensional fundamental solutions for transient heat transfer by conduction in an unbounded medium, half-space, slab and layered media, Eng Anal Bound Elem 30 (2006) 338–349. [86] 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 Bound Elem 45 (2014) 62–67. [87] H.W. Zhang, B.S. Chen, Y.X. Gu, An adaptive algorithm of precise integration for transient analysis, Acta Mech Solida Sin 14 (2001) 215–224. [88] C. Moler, C.V. Loan, Nineteen dubious ways to compute the exponential of a matrix, SIAM Rev 20 (1978) 801–836. [89] C. Moler, C.V. Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Rev 45 (2003) 1–46. [90] J.S. Agapiou, M.F. Devries, An experimental determination of the thermal conductivity of a 304L stainless steel powder metallurgy material, J Heat Tran 111 (1989) 281–286.

[61] K. Erharta, E. Divob, A.J. Kassaba, A parallel domain decomposition boundary element method approach for the solution of large-scale transient heat conduction problems, Eng Anal Bound Elem 30 (2006) 553–563. [62] W.Z. Feng, K. Yang, M. Cui, X.W. Gao, Analytically-integrated radial integration BEM for solving three-dimensional transient heat conduction problems, Int Commun Heat Mass Transfer 79 (2016) 21–30. [63] G.E. Hitti, M. Nemer, K.E. Khoury, Reducing CPU time for radiative exchange and transient heat transfer analysis using zone method, Numer Heat Transfer B 56 (2009) 23–37. [64] E.L. Wilson, R.E. Nickell, Application of the finite element method to heat conduction analysis, Nucl Eng Des 4 (1966) 276–286. [65] R. Cook, D. Malkus, M. Plesha, R. Witt, Concepts and applications of finite element analysis, J Pressure Vessel Technol 106 (2002) 188–194. [66] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, The finite element method: its basis and fundamentals, Butterworth-Heinemann, Oxford, 2005. [67] E.A. Thornton, Thermal structures for aerospace applications, VA, reston, (1996). [68] R.W. Lewis, K. Morgan, H. Thomas, K. Seetharamu, The finite element method in heat transfer analysis, John Wiley & Sons Inc., New York, 1996. [69] W. Wood, Practical time-stepping schemes, Claredon Press, Oxford, 1990. [70] W.Z. Feng, X.W. Gao, An interface integral equation method for solving transient heat conduction in multi-medium materials with variable thermal properties, Int J Heat Mass Transfer 98 (2016) 227–239. [71] F. Wu, Q. Gao, W.X. Zhong, Fast precise integration method for hyperbolic heat conduction problems, Appl Math Mech 34 (2013) 791–800. [72] Q. Gao, H.C. Cui, An efficient and accurate method for transient heat conduction in 1D periodic structures, Int J Heat Mass Transfer 108 (2017) 1535–1550. [73] B. Chen, Y. Gu, Z. Guan, H. Zhang, Nonlinear transient heat conduction analysis with precise time integration method, Numer Heat Transfer B 40 (2001) 325–341. [74] Q. Gao, W. Zhong, W.P. Howson, A precise method for solving wave propagation problems in layered anisotropic media, Wave Motion 40 (2004) 191–207. [75] Q. Gao, J.H. Lin, W.X. Zhong, F.W. Williams, Propagation of non-stationary random waves in viscoelastic stratified solids, Comput Geotech 33 (2006) 444–453. [76] Q. Gao, J.H. Lin, W.X. Zhong, W.P. Howson, F.W. Williams, Random wave propagation in a viscoelastic layered half space, Int J Solids Struct 43 (2006) 6453–6471. [77] J.H. Lin, W.P. Shen, F.W. Williams, Accurate high-speed computation of non-stationary random structural response, Eng Struct 19 (1997) 586–593.

231