A Novel Localized Meshless Method for Solving Transient Heat Conduction Problems in Complicated Domains

1  Introduction

Transient heat conduction exists widely in the thermal structure design of industrial equipment such as energy and chemical industry, as well as in the casting and heat treatment process of many new materials [1,2]. Therefore, the study on transient heat conduction is very important for a long time. It is not easy work to find analytical solutions to practical engineering problems. Hence, the numerical simulation becomes an important way of approximating numerical solutions.

There are many kinds of conventional numerical methods applied to transient heat conduction problems, for instance, the finite element method (FEM) [3,4], the finite difference method (FDM) [5] and the boundary element method (BEM) [6,7]. Among them, the FEM has been widely used in engineering and scientific calculation, and its mature commercial software has been widely used to analyze and simulate practical problems. There is no doubt that these traditional grid-type methods have their own advantages, but they still face many difficulties and challenges, such as complex mesh generation and other pre-processing processes, which affected the calculation accuracy and efficiency, especially for the complicated geometry.

In order to get rid of the complexity of mesh generation and reduce the time of preprocessing, various meshless methods have devoted considerable attention. These approaches include the element-free Galerkin method [8–11], the reproducing kernel particle method [12–15], the meshless local Petrov-Galerkin method [16,17], the radial basis function collocation method (RBFCM) [18,19], the generalized finite difference method (GFDM) [20–23], the singular boundary method (SBM) [24,25], the method of fundamental solutions (MFS) [26,27] and the boundary knot method (BKM) [28,29], etc. The successful application of these meshless methods fully demonstrates their development prospect. However, each method is faced with more or fewer limitations while playing its own advantages. As the semi-analytical, strong-form, and boundary-type meshless collocation methods, they are difficult to apply to complex-geometry and large-scale problems due to their dense matrix characteristics.

For the above reasons, a new class of meshless collocation techniques, called the localized semi-analytical meshless methods, has been proposed to solve various mathematical and mechanical problems. Such methods mainly include the localized method of fundamental solution (LMFS) [30], the localized singular boundary method (LSBM) [31], the localized Trefftz method (LTM) [32], and the localized boundary knot method (LBKM) [33]. The LMFS and the LSBM use the singular fundamental solutions as the basis functions, the LTM employs the T-complete functions, and the LBKM uses the non-singular general solutions. The first two methods need to determine the fictitious boundary and the source intensity factor resulting from the source singularity, while the latter two methods can be used for numerical approximation directly. Recently, these schemes have been successfully applied to heat and mass transfer [34,35], acoustics [36], elastic mechanics [37,38], inverse problems [39,40] and other aspects. For details, see [41] and references therein.

Very recently, an improved LBKM called the local knot method (LKM) [42–44], is developed to solve the convection-diffusion-reaction and acoustic problems. The improved scheme is a local RBF method, and does not need the boundary of supporting subdomain and artificial nodes on it, which enhances the computational efficiency. Compared with other local semi-analytical meshless methods, this approach adopts the non-singular general solution instead of singular fundamental solution, and thus it avoids the evaluation of fictitious boundary and the original intensity factor. Furthermore, the heat transfer equation becomes a modified Helmholtz equation after time approximation, whose solution can be easily solved by using the LKM. Considering the simplicity, accuracy, efficiency, and large-scale simulation of the LKM, this study employed this approach to solve transient heat conduction problems.

Firstly, the temporal derivative is discretized by using the finite difference scheme. The governing equation is transformed into a non-homogeneous modified Helmholtz equation. Secondly, the solution of transformed equation is decomposed into a particular solution and a homogeneous solution of the corresponding homogeneous equation. Finally, the particular solution and the homogeneous solution are obtained by the dual reciprocity method (DRM) [28,45] and the LKM, respectively. The particular solution only satisfies the non-homogeneous equation, rather than boundary condition. However, the homogeneous solution satisfies both homogeneous equation and corresponding boundary conditions. It should be noted that that there are many ways to approximate a particular solution of non-homogeneous modified Helmholtz equation besides the DRM, such as the multiple reciprocity method, the polynomial approximation, and the radial basis function approximation. Among them, the DRM is a simple and popular technique addressing non-homogeneous terms and therefore is used in this study.

The outline of this article is as follows. Section 2 first briefly reviews the basic theory of the transient heat conduction problem, and then introduces the LKM and DRM. In Section 3, three numerical examples are presented to demonstrate the accuracy and validity of the proposed scheme. Finally, Section 4 provides some conclusions.

2  Mathematical Formulations

2.1 Governing Equation

A bounded domain Ω with boundary Γ=∂Ω is considered. The governing equation of the transient heat conduction problems is indicated as

k∇2T(x,t)+Q(x,t)=ρc∂T(x,t)∂t,x∈Ω⊂Rd(1)

with initial and boundary conditions

T(x,t)=T¯(x,t),x∈ΓD(2)

−k∇T(x,t)⋅n=q¯(x,t),x∈ΓN(3)

T(x,0)=T0,x∈Ω(4)

where d denotes the dimension of space, T(x,t) the temperature on the point x at the time t, k the thermal conductivity, ρ the mass density, c the specific heat capacity, n the outward unit normal vector, Q the internal heat source, T0 the initial temperature. T¯(x,t) and q¯(x,t) represent the known temperature and heat flux on the boundaries ΓD and ΓN, respectively.

2.2 Time-Advancing Scheme

The transient heat conduction problem is transformed into the modified Helmholtz equation by using the finite difference scheme to deal with the time variable. For the time span [tn,tn+1]⊂[0,t]. The temperature T(x,t) at any node in the object, its partial derivative with respect to time and the internal heat source can be approximated by

T(x,t)=T(x,tn+1)(5)

∂T(x,t)∂t=T(x,tn+1)−T(x,tn)Δt(6)

Q(x,t)=Q(x,tn+1)(7)

where n denotes the time level, Δt=tn+1−tn is the time step.

Substituting Eqs. (5)–(7) into Eq. (1), the temperature Tn+1(x) at time tn+1 can be obtained. The sorted formula is as follows:

(∇2−ρckΔt)Tn+1(x)=−ρckΔtTn(x)−1kQn+1(x)(8)

It should be pointed out that the numerical solution Tn(x) in Eq. (8) can be obtained at time tn.

For convenience, Eq. (8) can be simplified as

(∇2−λ2)T(x)=f(x)(9)

in which

λ=ρckΔt(10)

f(x)=−ρckΔtTn(x)−1kQn+1(x)(11)

Noted that Eq. (9) is a non-homogeneous modified Helmholtz equation. Its solution can be decomposed into two parts (a particular solution and a homogeneous solution), namely,

T(x)=Tp(x)+Th(x)(12)

The particular solution only satisfies the non-homogeneous equation, rather than boundary condition. However, the homogeneous solution satisfies both homogeneous equation and corresponding boundary conditions. The next two subsections will introduce the dual reciprocity method and the LKM for approximating particular solution and homogeneous solution, respectively.

2.3 Dual Reciprocity Method

Obviously, the particular solution Tp(x) should satisfy non-homogeneous equation Eq. (9), namely,

(∇2−λ2)Tp(x)=f(x),x∈Ω(13)

This study uses the dual reciprocity method to approximate the particular solution Tp(x). According to its basic idea, the right-hand term of Eq. (13) is expressed by [28]

f(x)=∑j=1nαjφ(||x−yj||),x,yj∈Ω(14)

in which {αj}j=1n are undetermined coefficients, φ(⋅) is the radial basis function, {yj} is the interpolation points in the calculation, and n indicates the number of interpolation points.

Accordingly, the particular solution Tp(x) can be given by

Tp(x)=∑j=1nαjψ(||x−yj||)(15)

where ψ satisfies

(∇2−λ2)ψ(||x−yj||)=φ(||x−yj||)(16)

As we all know, there are many kinds of radial basis functions. The present study adopts the following function [28]:

φ(||x−yj||)=rj(17)

Substituting Eq. (17) into Eq. (16), we can get the expression of ψ

{ψ(||x−yj||)=−2λ4rj−rjλ2+2e−λrjλ4rj,rj>0ψ(||x−yj||)=−2λ3,rj=0(18)

where rj=||x−yj|| denotes the Euclidean distance between interpolation nodes.

2.4 Local Knot Method (LKM)

In this subsection, the LKM is used to approximate the homogeneous solution Th(x). First, we give the corresponding governing equation and boundary conditions

(∇2−λ2)Th(x)=0,x∈Ω(19)

Th(x)=T¯(x)−Tp(x),x∈ΓD(20)

k∂Th(x)∂n=q¯(x)−k∂Tp(x)∂n,x∈ΓN(21)

In the LKM, N=ni+nb1+nb2 nodes are discretized over the computational domain Ω, where ni, nb1 and nb2 are the numbers of interior nodes, Dirichlet boundary nodes and Neumann boundary nodes, respectively. Fig. 1 shows the schematic diagram of the LKM. Considering an arbitrary node (here called the central node x(0)) distributed in the computational domain, its m supporting nodes x(i),i=1,2,…,m can be easily determined. These nodes are contained in a circular or spherical subdomain called supporting domain. With regard to the m+1 local nodes inside the supporting subdomain Ωs, we can express the homogeneous solution as

Th(x(i))=∑j=0mςjG(rij),x(i)∈Ωs,i=0,1,…,m(22)

or for brevity

T(i)=G(i)ς(23)

where rij=||x(i)−x(j)||2, (ςj)j=0m are the unknown coefficients, G(rij) is the nonsingular general solution that can be expressed as [33]

G(rij)={I0(λrij),2Dsinh⁡(λrij)rij,3D(24)

where I0 is zero-order modified Bessel function of the first kind.

Figure 1: Schematic diagram of the LKM: (a) 2D; (b) 3D

Using the moving least squares theory, we define the following residual function:

ϕ(T)=∑i=0m[G(i)ς−T(i)]2ω(i)(25)

where ω(i) denotes the weighting function. In our calculation, the spline weight function [16] is used

ω(i)=1−6(didmax)2+8(didmax)3−3(didmax)4,i=0,1,…,m(26)

where di=||x(0)−x(i)||2, dmax=maxi=1,2,…,m(di) is the radius of supporting domain Ωs. According to Eq. (26), we can observe the smaller distance di make a much greater contribution for the LKM approximation.

The unknown coefficients ς=(ς0,ς1,…,ςm)T can be obtained by minimizing the function ϕ(T) with regard to ς, namely,

∂ϕ(T)∂ςj=0,j=0,1,…,m(27)

Then, a matrix equation can be obtained,

Aς=b(28)

where

A=[∑i=0mGi02ω(i)∑i=0mGi0Gi1ω(i)∑i=0mGi0Gi2ω(i)⋯∑i=0mGi0Gimω(i)∑i=0mGi12ω(i)∑i=0mGi1Gi2ω(i)⋯∑i=0mGi1Gimω(i)∑i=0mGi22ω(i)⋯∑i=0mGi2Gimω(i) SYM ⋱⋮∑i=0mGim2ω(i)],b=[∑i=0mGi0ω(i)T(i)∑i=0mGi1ω(i)T(i)∑i=0mGi2ω(i)T(i)⋮∑i=0mGimω(i)T(i)](29)

Further, the vector b in Eq. (29) can be recast as

b=[G00ω(0)G10ω(1)⋯Gm0ω(m)G01ω(0)G11ω(1)⋯Gm1ω(m)⋮⋮⋱⋮G0mω(0)G1mω(1)⋯Gmmω(m)](T(0)T(1)⋮T(m))=B(T(0)T(1)⋮T(m))(30)

In terms of Eqs. (28)–(30), ς=(ς0,ς1,…,ςm)T can be obtained by

ς=[ς0⋮ςm]=A−1B(T(0)T(1)⋮T(m))(31)

Let i=0 and using Eq. (23), we have

T(0)=G(0)ς=G(0)A−1B(T(0)T(1)⋮T(m))=∑j=0mα(j)T(j)(32)

or

T(0)−∑j=0mα(j)T(j)=0(33)

in which (α(j))j=0m=G(0)A−1B. Accordingly, the temperature at each node in the computational domain can be expressed in the same form as Eq. (33).

In addition, the normal derivative can be calculated by

∂T(0)∂nx(0)=∑j=0mςj∂G(x(0),x(j))∂nx(0)=∂G(0)∂nx(0)ς=∂G(0)∂nx(0)A−1B(T(0)T(1)⋮T(m))(34)

or for brevity

∂T(0)∂nx(0)=∑l=0mβ(l)T(l)(35)

where (β(l))l=0m=∂G(0)∂nx(0)A−1B.

At internal nodes and boundary nodes, the following equations should be satisfied:

T(i)−∑j=0mαi(j)T(j)=0,i=1,2,…,ni(36)

T(i)=T¯(i)−Tp(x(i)),i=ni+1,ni+2,…,ni+nb1(37)

∑l=0mβi(l)T(l)=q¯(x(i))−∂Tp(x(i))∂nx(i),i=ni+nb1+1,ni+nb1+2,…,N(38)

By combining Eqs. (36)–(38), we can obtain a linear system

ETh=h(39)

where EN×N is a sparse coefficient matrix, Th=[T(1),T(2),…,T(N)]T contains the homogeneous solutions at all nodes, and hN×1 is a vector consisting of right-hand term of and boundary data. After solving the above equation, the temperature T(x) at all nodes can thereby be obtained by the following formula:

T(x)=Tp(x)+Th(x)=∑j=1nαjψ(||x−yj||)+Th(x),x∈Ω(40)

3  Numerical Results and Discussion

Three numerical examples are investigated to test the accuracy and effectiveness of the proposed scheme. To measure numerical accuracy, the following errors are adopted:

Globalerror=∑k=1W(unk−uek)2/∑k=1W(uek)2(41)

Maxerror=max{|unk−uek|}(42)

in which W is the total number of test points, unk(k=1,2,…,W) and uek(k=1,2,…,W) are numerical and exact solutions at the test point, respectively. It should be pointed out that the nodes in examples 1 and 2 are generated from MATLAB software and the nodes in example 3 are automatically generated from Hypermesh software.

3.1 Example 1: Rectangular Domain with a Hole

In the first example, we consider a transient heat transfer problem in a rectangular plate with a hole [46]. The size parameters and boundary conditions of this model are given in Fig. 2. The density, the specific heat and the heat conductivity are ρ=5000 kg/m3, c=200J/(kg⋅°C) and k=200W/(m⋅°C), respectively. The initial temperature is taken to be T0=20°C, the temperature on the left boundary rises to T=500°C, while the right boundary maintains the initial temperature. Other boundaries are adiabatic. The temperatures at points A (0.07, 0.025) and B (0.075, 0.02) need to be studied and discussed in this example. Fig. 3 shows the distribution of nodes in the computation. The interior of plate is discretized into 3853 nodes and its boundary is discretized into 200 outer boundary nodes and 31 inner boundary nodes.

Figure 2: Dimensions and boundary conditions of the plate

Figure 3: Nodal distribution of the LKM

In the calculation, we set m=30 and Δt=4s. Since no analytical solution is available in this example, the finite element results calculated by COMSOL Multiphysics 5.4 are taken as the reference solution, in which the model is divided into 3853 domain elements. Fig. 4 gives the comparisons of the calculated temperatures between the FEM and the LKM at the points A and B within 200 s. From Fig. 4, it can be seen that the LKM results are in good agreement with those of the FEM, and the temperature at point A is higher than that at point B. In addition, when the time is less than 10 s, there is a certain deviation between the results of the LKM and the FEM, but with the increase of time, the deviation between them decreases and the temperature tends to be stable.

Figure 4: Temperature variation at points A and B with time

Then, we consider the numbers of supporting points as m=20,30,40,50 to verify the accuracy and convergence of the LKM. Fig. 5 shows the relative deviation curves of the FEM and the proposed method at point A. It can be clearly seen from Fig. 5 that the deviations decrease rapidly at first, reach the optimal value and then tend to be stable as time goes on. When the number of supporting points is 20 and the time is 100 s, the minimum deviation is 1.148×10−6. At the 200 s, the deviations of the four cases are all less than 4.769×10−3, meeting the accuracy requirements. Furthermore, with the increase of supporting points, the relative deviations converge faster and reach the optimal value faster, indicating that LKM has good convergence.

Figure 5: The relative deviation curves at point A

Figs. 6–8 illustrate the temperature distributions calculated by the LKM and the FEM at 4, 8 and 200 s, respectively. Obviously, the results near the hole have a little different at the beginning of the time period, but gradually become identical when the time goes to 200 s.

Figure 6: Profiles of temperature on the plate at time 4 s: (a) FEM; (b) LKM

Figure 7: Profiles of temperature on the plate at time 8 s: (a) FEM; (b) LKM

Figure 8: Profiles of temperature on the plate at time 200 s: (a) FEM; (b) LKM

3.2 Example 2: A Cube with Heat Generation

The LKM is applied to analyze the temperature field of the cubic domain in the second example. The shape parameters and initial temperature of the cube are shown in Fig. 9a. Fig. 9b shows the distribution of nodes in the computation, including 5832 boundary nodes and 1944 internal nodes. The thermal diffusivity is assumed to be a=k/ρc=0.16m2/s. The heat generation is expressed as [28]

Q(x,y,z)=Q0sin⁡πxLsin⁡πyLsin⁡πzL(t≥0)(43)

Figure 9: Dimensions and nodal distribution of the cube: (a) Dimensions of the cube; (b) Distribution of nodes

The analytical solution is given by

T(x,y,z)=L2Q03π2ksin⁡πxLsin⁡πyLsin⁡πzL[1−exp⁡(−3aπ2tL2)](44)

where Q0/k=100°C/m2.

In this example, we consider two cases: (a) All the boundaries satisfy the Dirichlet boundary conditions; (b) The boundary y=0 satisfies the Neumann boundary conditions, and the heat flux can be expressed as q¯(x,t)=−LQ03πsin⁡πxLsin⁡πzL[1−exp⁡(−3aπ2tL2)], while the temperature is 0°C on the remaining boundaries.

In this example, the number of the supporting nodes is 60, and the time interval is Δt=0.05s. For the two different cases, Fig. 10 shows the comparison of temperature on the line {(x,y,z)|0≤x≤1,y=z=0.5} at different times obtained by using the LKM and analytical expression. It can be plainly indicated that the developed method has good accuracy for the simulation of transient heat conductions with two kinds of boundary condition. In addition, it can also be noted that the error is slightly large when x∈[0.4 0.6], but the error gradually decreases along with the time goes on.

Figure 10: Temperature distributions on the line {(x,y,z)|0≤x≤1,y=z=0.5}: (a) Case 1;  (b) Case 2

Additionally, Fig. 11 plots the temperature curves with time marching for three internal points C (0.2, 0.2, 0.2), D (0.4, 0.4, 0.4) and E (0.6, 0.6, 0.6) for two cases and compares with the corresponding exact results. The results show that the method has high precision.

Figure 11: Distributions of temperature at different times at points C, D and E: (a) Case 1; (b) Case 2

By comparing the results of the LKM and the BKM in different cases, the accuracy and convergence of this method are carefully validated. With the increase of the total number of nodes, the variation of the errors and the condition numbers of the two methods are listed in Tables 1 and 2. As can be seen from Tables 1 and 2, although there are a few errors at some nodes, the results are still clear that the maximum and global errors of the two methods decrease with the increase of the number of nodes, and meet certain convergence trend. Next, we can conclude that when the total number of nodes is greater than 3211, the LKM is better than the BKM. Moreover, the LKM has fewer condition numbers regardless of the total number of nodes. In summary, the proposed LKM is high-precision, stable and convergent for solving 3D transient heat transfer problems.

Finally, we can observe from Fig. 12 that the distributions of absolute errors on the cross section z=0.5 in two cases are different. In case 1, it can be clearly seen that the absolute error in the center of the computational domain is largest (but less than 1.711 × 10−2), and gradually decreases from the center to the boundaries. In case 2, we can find that the absolute error on the boundary y=0 is largest (but less than 7.21 × 10−2), and decreases sharply from it to other boundaries. The above phenomenon is mainly caused by different boundary conditions.

Figure 12: Distributions of absolute error on the cross section (z = 0.5) for (a) Case 1; (b) Case 2

3.3 Example 3: A Hollow Cylinder

In order to verify the proposed LKM for simulating 3D transient heat conduction in a multi-connected domain, the third example considers is a hollow cylinder. Fig. 13a plots the geometry parameters and boundary conditions of the model. The initial temperature is set to be T0=500°C. The temperature at the top surface is suddenly decreased to 50°C, while the temperature stays the same at the bottom surface, and other surfaces are thermally insulated. The density, the heat conductivity and the specific heat are ρ=7850 kg/m3, k=60W/(m⋅°C) and c=460J/(kg⋅°C), respectively. The discretization of nodes in the computation is shown in Fig. 13b. There are 14146 nodes in the computation, of which 9706 are boundary nodes and the rest are internal nodes. Then, we continue to set FEM solution with more nodes (includes 229323 domain elements) than the proposed method as the reference solution.

Figure 13: Dimensions, boundary conditions and nodal distribution of the hollow cylinder: (a) Dimensions and boundary conditions, (b) Nodal distribution

In this case, the number of supporting nodes is set to be 60. The length of each time step is Δt=1s. Figs. 14 and 15 depict the distribution of temperature field on the surface of the computational domain at the time 10 and 100 s. From Figs. 14 and 15, it can be observed that the temperature increases gradually from top to bottom of the hollow cylinder. When the time advances to 100 s, the temperature changes in the calculation domain tends to be stable, and the results are almost identical for both the LKM and the FEM. We can carefully conclude that the LKM can also have high accuracy when using fewer nodes compared with the FEM. A large number of numerical experiments show that the proposed method has great potential in solving transient heat transfer problems, but the method needs to be further improved through continuous theoretical verification and engineering application.

Figure 14: Contours of temperature at time 10 s: (a) LKM; (b) FEM

Figure 15: Contours of temperature at time 100 s: (a) LKM; (b) FEM

4  Conclusion

This paper presented a novel coupling algorithm of the LKM and the dual reciprocity method to simulate the transient heat transfer numerically. The LKM is a semi-analytical and local meshless method which uses the non-singular general solution as the basis function of interpolation. Unlike the existing conventional methods, the coupled approach avoids mesh division, reduces condition number and expands the simulation scale.

Numerical results of the three benchmark examples involving simply-connected and multiply-connected domains have been investigated in detail. In the examples without exact solutions, our numerical results are compared with the FEM results from COMSOL software. Numerical experiments indicate that the proposed scheme is accurate and effective, and has small condition numbers.

This paper focuses on the application of the LKM in transient heat transfer problems, and shows its superior performance in solving such problems. However, there are still some issues that need to be addressed. The present study uses the Euler formula to discretize the time derivative, which has a certain impact on the accuracy and efficiency of the present method. Some more complex and accurate methods addressing the time derivative are interesting and will be discussed in more detail subsequently. Furthermore, the proposed methodology can be further optimized and improved to accurately solve the high-dimensional nonlinear heat transfer as well as the coupling heat transfer in complicated geometry.

Funding Statement: This work is supported by the National Natural Science Foundation of China (No. 11802151), the Natural Science Foundation of Shandong Province of China (No. ZR2019BA008), and the China Postdoctoral Science Foundation (No. 2019M652315).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.