A Novel Hybrid Boundary-Type Meshless Method for Solving Heat Conduction Problems in Layered Materials

In this article, we propose a novel meshless method for solving two-dimensional stationary heat conduction problems in layered materials. The proposed method is a recently developed boundary-type meshless method which combines the collocation scheme from the method of fundamental solutions (MFS) with the collocation Trefftz method (CTM) to improve the applicability of the method for solving boundary value problems. Particular non-singular basis functions from cylindrical harmonics are adopted in which the numerical approximation is based on the superposition principle using the non-singular basis functions expressed in terms of many source points. For the modeling of multi-layer composite materials, we adopted the domain decomposition method (DDM), which splits the domain into smaller subdomains. The continuity of the flux and the temperature has to be satisfied at the interface of subdomains for the problem. The validity of the proposed method is investigated for several test problems. Numerical applications were also carried out. Comparison of the proposed method with other meshless methods showed that it is highly accurate and computationally efficient for modeling heat conduction problems, especially in heterogeneous multi-layer composite materials.


Introduction
Since the 1940s, the knowledge of using composite laminate materials in industry has improved significantly [1].A variety of methods have been developed to deal with composite laminate materials, which make the subject of composite laminate materials a matured discipline in applied science at present.The industrial use of composite laminate materials is widespread because they have the advantages of their anisotropic nature, which allows the material to be applied to a variety of engineering applications.However, composite materials may fail if subjected to severe environments such as high temperatures, even though there is no external load applied to the composites [2].Specifically, failure in layered composites is caused by thermal induced stresses which are generated at the interface between different materials because the temperature distribution in the composites is non-uniform or discontinuous [3].In order to understand how thermal stress is generated and distributed in the layered composites, it is important to investigate the temperature distribution in the composite laminate materials.While the mechanical behavior of composite laminate materials has been explored for decades, studies on the heat conduction problems of layered composites using the meshless methods are less explored, which initiated this study.
In the past, many heat conduction problems in layered composites have been routinely solved by numerical methods such as the finite difference method [4][5][6], finite element method [7,8], finite volume method [9][10][11] and the boundary element method [12,13].In contrast to mesh-based numerical methods, the meshless methods, which do not need the mesh generation and boundary integral have been proposed such as the analytical method [14], method of fundamental solutions (MFS) [15][16][17], boundary knot method (BKM) [18], collocation Trefftz method (CTM) [19][20][21], radial basis function collocation method (RBFCM) [22][23][24][25], element-free Galerkin method (EFG) [26], reproducing kernel particle method (RKPM) [27,28], modified polynomial expansion method [29], meshless local boundary integral equation method (LBIE) [30,31], and so on.Among these, boundary-type meshless methods have attracted considerable attention because of their simplicity.The CTM [32] is categorized as a boundary-type meshless method for the numerical solution of problems where approximate solutions are expressed as a truncated series of T-complete basis functions automatically satisfying governing equations.The use of the CTM is less widespread because the system of linear equations obtained from the Trefftz method is an ill-posed system [33].On the other hand, the MFS is also a boundary-type meshless method for solving problems where the solutions are approximated by the fundamental solution which is expressed in terms of source points.The MFS requires placing source points outside the domain of the problem to avoid the effects of the singular characteristics of the fundamental solution; however, it often encounters difficulties such as finding an appropriate location for the source points [34][35][36].
In this study, we propose a novel hybrid boundary-type meshless method for solving two-dimensional stationary heat conduction problems in layered composite materials.The proposed method combines the collocation scheme from the MFS with the CTM to improve the applicability of both methods.Particular non-singular basis functions from the cylindrical harmonics are adopted in which the numerical solutions are approximated by superpositioning of the non-singular basis functions expressed in terms of many source points.For the modeling of multi-layer composite materials, we adopted the domain decomposition method (DDM) [37], which splits the domain into smaller subdomains that are intersected only at the interface between layers.For each subdomain, there exists an independent layer with its own thermal conductivity.The basic concept of the DDM for modeling composite materials is that the continuity of the flux and the temperature has to be satisfied at the interface of subdomains for the problem.The validity of the proposed method is investigated for several test problems.Numerical analysis is also carried out.The rest of this article is organized as follows.In Section 2, we describe the governing equation of two-dimensional stationary heat conduction problems.Section 3 is devoted to the formulation of the hybrid boundary-type meshless method.In Section 4, numerical analysis of several test problems were conducted to evaluate the performance of the proposed numerical scheme.Finally, our conclusions are presented in Section 5.

Problem Statement
Considering a two-dimensional domain, Ω, enclosed by a boundary Γ, the governing equations of stationary heat conduction and boundary conditions can be expressed as follows: and u = g on Γ D (2) where u denotes the temperature, ∇ denotes the Laplace operator, Ω denotes the object domain under consideration, n denotes the outward normal direction, g denotes the boundary where the Dirichlet

Hybrid Boundary-Type Meshless Method for Modeling Heat Conduction Problems
In the polar coordinate system, the governing equation of the stationary heat conduction problems is shown as follows: where ρ is the radius and θ is the polar angle in the polar coordinate system.In the proposed method, the approximation is established via superpositioning a series of particular non-singular basis functions of the heat equation in terms of many source points located within the domain.
An approximation solution of two-dimensional stationary heat conduction problems can be expressed as a linear combination of particular non-singular basis function with unknown coefficients, a j,k .
where x = (ρ, θ) and x ∈ ∂Ω is the spatial coordinate that is collocated on the boundary, y j = (ρ j , θ j ) is the source point, O is the number of source points, P is the order of the particular non-singular basis function, a j,k = [ c jk d jk ] is a vector of unknown coefficients to be determined, H(x, y j ) is the selected combination of the particular solutions of the heat equation in the two-dimensional polar coordinate system.The particular non-singular basis function of the heat equation can be expressed as In this study, we adopt the characteristic length [38] as where max(ρ) is the maximum radial distance in the problem domain.Applying the Dirichlet boundary condition, we obtain where g(x l ) is the Dirichlet boundary value imposed at the boundary collocation point.Q is the number of boundary collocation points.The Neumann boundary condition can be written as where where n x and n y are normal vectors in the x and y direction, respectively.∂u/∂ρ and ∂u/∂θ can be rewritten as follows: Inserting Equations ( 11) and ( 12) into Equation ( 9), we may obtain where Applying the Neumann boundary condition, we obtain where f (x l ) is the Neumann boundary value imposed at boundary collocation points.Using Equations ( 8) and ( 16), we obtain the following linear systems of the form: ... where where A is a matrix of particular non-singular basis function with the size of Q × S, S = (2P × O), Q is the number of boundary collocation points for given boundary conditions, P is the order of the particular non-singular basis function, and O is the number of source points.α is a vector of unknown coefficients with the size of S × For the modeling of multi-layer composite materials, we adopted the DDM [37].The DDM splits the domain into smaller subdomains in which subdomains are intersected only at the interface between layers.For each subdomain, there exists an independent layer with its own thermal conductivity.The boundary and source points are collocated in each subdomain.At the interface, the boundary collocation points on the left and right sides coincide with each other.The basic concept of the DDM for modeling composite materials is that the subdomains intersecting at the interface must satisfy flux conservation and continuity of the temperature for the problem, so that flux and temperature at the interface between two consecutive layers remain the same.For instance, we consider a rectangular domain Ω, which can be divided into two smaller subdomains, Ω 1 and Ω 2 as shown in Figure 1.In order to simulate the heat conduction problem, the rectangular domain is divided into Γ 1 , Γ 2 , . . ., Γ 8 sub-boundaries.At Ω 1 subdomains, the sub-boundaries include Γ 1 , Γ 2 , Γ 3 and Γ 4 ; At Ω 2 subdomains, the sub-boundaries include Γ 5 , Γ 6 , Γ 7 and Γ 8 .As noted above, the sub-boundary at the interface should satisfy the flux conservation and continuity of the temperature between two consecutive materials.The additional boundary conditions at the interface can be expressed as Matching the Dirichlet and Neumann boundary conditions on boundary collocation points, we may obtain the following simultaneous linear equations where A Ω 1 with the size of l 1 × P 1 and A Ω 2 with the size of l 2 × P 2 are the A matrix shown in Equation (18) for Ω 1 and Ω 2 , respectively.l 1 and l 2 are the number of boundary collocation points for Ω 1 and Ω 2 , respectively.P 1 and P 2 are the number of the terms related to the order of the particular non-singular basis function for Ω 1 and Ω 2 , respectively.A I | Γ 2 of the boundary Γ 2 with the size of l I × P 1 and A I | Γ 6 of the boundary Γ 6 with the size of l I × P 2 are the A matrices at the interface.l I is the number of boundary collocation points at the interface.0 Ω 1 and 0 Ω 2 are zero matrices with the size of l 2 × P 1 and l 1 × P 2 , respectively.α Ω 1 with the size of P 1 × 1 is a vector of unknown coefficient of Ω 1 , α Ω 2 with the size of P 2 × 1 is a vector of an unknown coefficient of Ω 2 .α Ω 1 and α Ω 2 are vectors of boundary values of Ω 1 and Ω 2 , respectively.b I = [0 g 0 f ] T , 0 g and 0 f are zero vectors with the size of l I × 1 for Dirichlet and Neumann boundary conditions at the interface, respectively.By solving the above simultaneous linear equations, we may obtain two sets of unknown coefficients, b Ω 1 and b Ω 2 , for subdomains, Ω 1 and Ω 2 , respectively.To obtain the temperature for subdomains, the inner collocation points in Ω 1 and Ω 2 must be placed.The temperature, u, at inner collocation points can then be found by Equation (5) using A Ω 1 and b Ω 1 for Ω 1 , as well as A Ω 2 and b Ω 2 for Ω 2 .

Numerical Examples
To validate the proposed hybrid boundary-type meshless method for modeling stationary heat conduction problems, several test problems were investigated during the simulation of heat conduction problems in composite materials and compared with the analytical solutions.In order to evaluate the accuracy of the numerical solution, we used the absolute and relative error as shown in the following equations.
)) ( ( where are the exact and numerical solutions for the inner points i x , respectively.

Example 1
The first scenario under investigation is a two-dimensional homogeneous isotropic stationary heat conduction problem [39], as shown in Figure 2.For a two-dimensional domain, Ω , enclosed by an elliptical-shaped section, the governing equation of the heat conduction problem is expressed as follows: The two-dimensional object boundary under consideration is defined as The profile of the composite materials for the analysis.

Numerical Examples
To validate the proposed hybrid boundary-type meshless method for modeling stationary heat conduction problems, several test problems were investigated during the simulation of heat conduction problems in composite materials and compared with the analytical solutions.In order to evaluate the accuracy of the numerical solution, we used the absolute and relative error as shown in the following equations.
Absolute error = u exact (x i ) − u nu (x i ) Relative error = (u exact where u exact (x i ) and u nu (x i ) are the exact and numerical solutions for the inner points x i , respectively.

Example 1
The first scenario under investigation is a two-dimensional homogeneous isotropic stationary heat conduction problem [39], as shown in Figure 2.For a two-dimensional domain, Ω, enclosed by an elliptical-shaped section, the governing equation of the heat conduction problem is expressed as follows: The two-dimensional object boundary under consideration is defined as The analytical solution may be found as For solving this example, the order of the particular non-singular basis function P, is 15.The boundary condition is assumed to be the Dirichlet boundary condition using the analytical solution which is expressed as Equation ( 30) for the problem.The boundary and source points are uniformly distributed on the boundary for the hybrid boundary-type meshless method, as demonstrated in Figure 3.We compared the computed temperature of the inner collocation points in the computational domain with those from the analytical solution, as shown in Figure 4a.We obtained highly accurate numerical solutions in the order of 10 −16 for this problem.Figure 4b demonstrates the accuracy of the element-free Galerkin method (EFG) [39] and the improved interpolating element-free Galerkin (IIEFG) method [39].It is found that the accuracy of the EFG and the IIEFG methods can only reach the order of 10 −16 for the same problem.Figure 5 shows the comparison of the computed results with the analytical solution of temperature along line AB.In addition, Figure 6 depicts a plot of the comparison of the computed temperature distribution with the analytical solution.From the figures, it can be seen that the numerical solution agrees very well with the analytical solution.This example shows that the proposed method may obtain more accurate results than the EFG and IIEFG methods.
The analytical solution may be found as For solving this example, the order of the particular non-singular basis function P , is 15.The boundary condition is assumed to be the Dirichlet boundary condition using the analytical solution which is expressed as Equation ( 30) for the problem.The boundary and source points are uniformly distributed on the boundary for the hybrid boundary-type meshless method, as demonstrated in Figure 3.We compared the computed temperature of the inner collocation points in the computational domain with those from the analytical solution, as shown in Figure 4a.We obtained highly accurate numerical solutions in the order of 10 −16 for this problem.Figure 4b demonstrates the accuracy of the element-free Galerkin method (EFG) [39] and the improved interpolating element-free Galerkin (IIEFG) method [39].It is found that the accuracy of the EFG and the IIEFG methods can only reach the order of 10 −16 for the same problem.Figure 5 shows the comparison of the computed results with the analytical solution of temperature along line AB.In addition, Figure 6 depicts a plot of the comparison of the computed temperature distribution with the analytical solution.From the figures, it can be seen that the numerical solution agrees very well with the analytical solution.This example shows that the proposed method may obtain more accurate results than the EFG and IIEFG methods.
The analytical solution may be found as For solving this example, the order of the particular non-singular basis function P , is 15.The boundary condition is assumed to be the Dirichlet boundary condition using the analytical solution which is expressed as Equation ( 30) for the problem.The boundary and source points are uniformly distributed on the boundary for the hybrid boundary-type meshless method, as demonstrated in Figure 3.We compared the computed temperature of the inner collocation points in the computational domain with those from the analytical solution, as shown in Figure 4a.We obtained highly accurate numerical solutions in the order of 10 −16 for this problem.Figure 4b demonstrates the accuracy of the element-free Galerkin method (EFG) [39] and the improved interpolating element-free Galerkin (IIEFG) method [39].It is found that the accuracy of the EFG and the IIEFG methods can only reach the order of 10 −16 for the same problem.Figure 5 shows the comparison of the computed results with the analytical solution of temperature along line AB.In addition, Figure 6 depicts a plot of the comparison of the computed temperature distribution with the analytical solution.From the figures, it can be seen that the numerical solution agrees very well with the analytical solution.This example shows that the proposed method may obtain more accurate results than the EFG and IIEFG methods.

Example 2
In order to verify the proposed hybrid boundary-type meshless method, we considered a homogeneous and isotropic example with a complex boundary shape.As shown in Figure 7, this example is a two-dimensional steady state heat conduction problem [39].The governing equation of the example enclosed by Ω can be expressed as The two-dimensional object boundary under consideration is defined as: The analytical solution may be found as For solving this example, the order of the particular non-singular basis function P , is 15.The boundary condition is assumed to be the Dirichlet boundary condition using the analytical solution which is expressed as Equation (37).The boundary and source points are uniformly distributed on the boundary for the hybrid boundary-type meshless method, as demonstrated in Figure 8.
We first compared the computed temperature of the inner collocation points in the computational domain with those from the analytical solution, as shown in Figure 9a.We obtained highly accurate numerical solutions in the order of 10 −11 for this problem.Figure 9b demonstrates the accuracy of the element-free Galerkin method (EFG) [39] and the improved interpolating element-free Galerkin (IIEFG) method [39].It was found that the accuracy of the EFG and the IIEFG methods can only reach to the order of 10 −3 for the same problem.Figure 10 shows the comparison of the computed results with the analytical solution of temperature along line OA.In addition, Figure 11 depicts a plot of the comparison of the computed temperature distribution with the analytical solution.From these figures, it can be observed that the numerical solution agrees very well with the analytical solution.Again, this example shows that the proposed method may obtain more accurate results than the EFG and IIEFG methods.

Example 2
In order to verify the proposed hybrid boundary-type meshless method, we considered a homogeneous and isotropic example with a complex boundary shape.As shown in Figure 7, this example is a two-dimensional steady state heat conduction problem [39].The governing equation of the example enclosed by Ω can be expressed as The two-dimensional object boundary under consideration is defined as: The analytical solution may be found as For solving this example, the order of the particular non-singular basis function P, is 15.The boundary condition is assumed to be the Dirichlet boundary condition using the analytical solution which is expressed as Equation (37).The boundary and source points are uniformly distributed on the boundary for the hybrid boundary-type meshless method, as demonstrated in Figure 8.
We first compared the computed temperature of the inner collocation points in the computational domain with those from the analytical solution, as shown in Figure 9a.We obtained highly accurate numerical solutions in the order of 10 −11 for this problem.Figure 9b demonstrates the accuracy of the element-free Galerkin method (EFG) [39] and the improved interpolating element-free Galerkin (IIEFG) method [39].It was found that the accuracy of the EFG and the IIEFG methods can only reach to the order of 10 −3 for the same problem.Figure 10 shows the comparison of the computed results with the analytical solution of temperature along line OA.In addition, Figure 11 depicts a plot of the comparison of the computed temperature distribution with the analytical solution.From these figures, it can be observed that the numerical solution agrees very well with the analytical solution.Again, this example shows that the proposed method may obtain more accurate results than the EFG and IIEFG methods.

Example 3
The third example under investigation is the modeling of two-dimensional isotropic stationary heat conduction problems in composite materials with a rectangular domain [40], as shown in Figure 1.The thermal conductivities in layer 1 and layer 2 are , respectively.
For a two-dimensional domain, Ω , enclosed by a rectangular boundary, the governing equation is expressed as follows: where k denotes thermal conductivity.The two-dimensional object boundary under consideration is defined as For the heat conduction equation, the analytical solution can be obtained as k are thermal conductivities in layer 1 and layer 2, respectively.The boundary and source points are uniformly distributed on the boundary, as demonstrated in Figure 12.For modeling the heat conduction problem in composite materials with a rectangular domain, the DDM [37] is adopted so that the domain boundary can be divide into several subdomains.Figure 1 shows that the two-dimensional object boundary is divided into several sub-boundaries: 1

Example 3
The third example under investigation is the modeling of two-dimensional isotropic stationary heat conduction problems in composite materials with a rectangular domain [40], as shown in Figure 1.The thermal conductivities in layer 1 and layer 2 are 1 W/m K and 2 W/m K, respectively.For a two-dimensional domain, Ω, enclosed by a rectangular boundary, the governing equation is expressed as follows: where k denotes thermal conductivity.The two-dimensional object boundary under consideration is defined as For the heat conduction equation, the analytical solution can be obtained as where k 1 and k 2 are thermal conductivities in layer 1 and layer 2, respectively.The boundary and source points are uniformly distributed on the boundary, as demonstrated in Figure 12.For modeling the heat conduction problem in composite materials with a rectangular domain, the DDM [37] is adopted so that the domain boundary can be divide into several subdomains.Figure 1 shows that the two-dimensional object boundary is divided into several sub-boundaries: The Dirichlet boundary condition is imposed on the sub-boundaries using the analytical solution as shown in Equation (41) for the problem.At the interface, the sub-boundaries, including Γ 2 and Γ 6 , should satisfy the flux conservation and the continuity of temperature.Therefore, the boundary conditions are given as follows: In this example, we utilized the hybrid boundary-type meshless method for modeling heat conduction problems in composite materials with a rectangular domain.The reported CPU time is 1.67 seconds using the computer with Intel Core i7-7700 CPU at 3.60 GHz and the equipment is manufactured by the ASUSTeK Computer Inc., Taipei, Taiwan.To validate the accuracy of the proposed method, we collocated 1200 inner points inside the domain uniformly to obtain the solution of the temperature.Figure 13 demonstrates the comparison of the computed results with those from the analytical solution.It is found that the results agree very well with each other.
To evaluate the accuracy of the proposed method, the absolute error and relative error of the computed results with the analytical solution were calculated as shown in Figure 14.We obtained highly accurate numerical solutions in the order of 10 −11 for this problem.The singular boundary method [40] has been used to model the same problem.As shown in Table 1, it was found that the accuracy of the singular boundary method for this example only reached to the order of 10 −3 .
Appl.Sci.2018, 8, x FOR PEER REVIEW 13 of 24 In this example, we utilized the hybrid boundary-type meshless method for modeling heat conduction problems in composite materials with a rectangular domain.The reported CPU time is 1.67 seconds using the computer with Intel Core i7-7700 CPU at 3.60 GHz and the equipment is manufactured by the ASUSTeK Computer Inc., Taipei, Taiwan.To validate the accuracy of the proposed method, we collocated 1200 inner points inside the domain uniformly to obtain the solution of the temperature.Figure 13 demonstrates the comparison of the computed results with those from the analytical solution.It is found that the results agree very well with each other.
To evaluate the accuracy of the proposed method, the absolute error and relative error of the computed results with the analytical solution were calculated as shown in Figure 14.We obtained highly accurate numerical solutions in the order of 10 −11 for this problem.The singular boundary method [40] has been used to model the same problem.As shown in Table 1, it was found that the accuracy of the singular boundary method for this example only reached to the order of 10 −3 .

Example 4
The fourth example under investigation is the modeling of two-dimensional isotropic stationary heat conduction problems of composite materials in a square domain with a hole [40], as depicted in Figure 15.The thermal conductivities in layer 1 and layer 2 are , respectively.In the first layer, two squares with sides f 2 and 20 around the boundary domain, respectively and the second layer is composed of two squares with sides of 20 and 60 around the boundary domain.For a two-dimensional domain, Ω , enclosed by a boundary with a square domain, the governing equation is shown in Equation (38).The analytical solution can be obtained as where k and 2 k are thermal conductivities in layer 1 and layer 2, respectively.The boundary and source points are uniformly distributed on the boundary, as demonstrated in Figure 16.For modeling the heat conduction problem with composite materials, the DDM [37] was adopted so that the domain boundary can be divided into several subdomains.Figure 15 shows that the two-dimensional object boundary is divided into several sub-boundaries: Therefore, the boundary conditions are given as follows.

Example 4
The fourth example under investigation is the modeling of two-dimensional isotropic stationary heat conduction problems of composite materials in a square domain with a hole [40], as depicted in Figure 15.The thermal conductivities in layer 1 and layer 2 are 1 W/m K and 2 W/m K, respectively.In the first layer, two squares with sides f 2 and 20 around the boundary domain, respectively and the second layer is composed of two squares with sides of 20 and 60 around the boundary domain.For a two-dimensional domain, Ω, enclosed by a boundary with a square domain, the governing equation is shown in Equation (38).The analytical solution can be obtained as where k 1 and k 2 are thermal conductivities in layer 1 and layer 2, respectively.The boundary and source points are uniformly distributed on the boundary, as demonstrated in Figure 16.For modeling the heat conduction problem with composite materials, the DDM [37] was adopted so that the domain boundary can be divided into several subdomains.Figure 15 shows that the two-dimensional object boundary is divided into several sub-boundaries: Γ 1 , Γ 2 , Γ 3 , . . ., Γ 16 .The Dirichlet boundary condition is imposed on the sub-boundaries using the analytical solution for the problem, as shown in Equation (43).At the interface, the sub-boundaries, including Γ 5 , Γ 6 , Γ 7 , Γ 8 , Γ 9 , Γ 10 , Γ 11 and Γ 12 , should satisfy the flux conservation and the continuity of temperature.Therefore, the boundary conditions are given as follows.The reported CPU time was 18.00 s using the Intel Core i7-7700 CPU at 3.60 GHz.To examine the accuracy of the proposed method, we collocated 5445 inner points inside the domain uniformly   The reported CPU time was 18.00 s using the Intel Core i7-7700 CPU at 3.60 GHz.To examine the accuracy of the proposed method, we collocated 5445 inner points inside the domain uniformly The reported CPU time was 18.00 s using the Intel Core i7-7700 CPU at 3.60 GHz.To examine the accuracy of the proposed method, we collocated 5445 inner points inside the domain uniformly to obtain the solution of the temperature.The computed temperature distribution of the inner points is depicted in Figure 17.It was found that the computed results agree very well with those from the analytical solution.To evaluate the accuracy of the proposed method, the absolute error and relative error of the computed results with the analytical solution were calculated as shown in Figure 18.We obtained highly accurate numerical solutions in the order of 10 −7 and 10 −11 for the absolute error and relative error, respectively.The singular boundary method [40] was used to model the same problem.As shown in Table 2, it was found that the accuracy of the singular boundary method for this example can only reach to the order of 10 −3 .
Appl.Sci.2018, 8, x FOR PEER REVIEW 17 of 24 to obtain the solution of the temperature.The computed temperature distribution of the inner points is depicted in Figure 17.It was found that the computed results agree very well with those from the analytical solution.To evaluate the accuracy of the proposed method, the absolute error and relative error of the computed results with the analytical solution were calculated as shown in Figure 18.We obtained highly accurate numerical solutions in the order of 10 −7 and 10 −11 for the absolute error and relative error, respectively.The singular boundary method [40] was used to model the same problem.As shown in Table 2, it was found that the accuracy of the singular boundary method for this example can only reach to the order of 10 −3 .to obtain the solution of the temperature.The computed temperature distribution of the inner points is depicted in Figure 17.It was found that the computed results agree very well with those from the analytical solution.To evaluate the accuracy of the proposed method, the absolute error and relative error of the computed results with the analytical solution were calculated as shown in Figure 18.We obtained highly accurate numerical solutions in the order of 10 −7 and 10 −11 for the absolute error and relative error, respectively.The singular boundary method [40] was used to model the same problem.As shown in Table 2, it was found that the accuracy of the singular boundary method for this example can only reach to the order of 10 −3 .

Example 5
In order to verify the hybrid boundary-type meshless method, we considered a more complex boundary shape with an irregular domain by comparing the results with the analytical solution.The fifth scenario investigated is similar to the fourth scenario, but it involved a more complex boundary shape, as depicted in Figure 19.The thermal conductivities in layer 1 and layer 2 are , respectively.The governing equation of Example 5 is expressed in Equation (35).
Considering a boundary with an irregular shape, we may write the equations of the two-dimensional boundary as follows: For the heat conduction equation, the analytical solution can be obtained as

Example 5
In order to verify the hybrid boundary-type meshless method, we considered a more complex boundary shape with an irregular domain by comparing the results with the analytical solution.The fifth scenario investigated is similar to the fourth scenario, but it involved a more complex boundary shape, as depicted in Figure 19.The thermal conductivities in layer 1 and layer 2 are 1 W/m K and 2 W/m K, respectively.The governing equation of Example 5 is expressed in Equation (35).Considering a boundary with an irregular shape, we may write the equations of the two-dimensional boundary as follows: Γ 10 = { (x 10 , y 10 )|x 10 = ρ(θ 10 ) cos θ 10 , y 10 = ρ(θ 10 ) sin θ 10 }, ρ(θ 10 ) = 20(e (sin θ 10 sin 2θ 10 )2 + e (cos θ 10 cos 2θ 10 )2 ), 0 ≤ θ 10 ≤ 2π (51) For the heat conduction equation, the analytical solution can be obtained as where k 1 and k 2 are thermal conductivities in layer 1 and layer 2, respectively.The boundary and source points are uniformly distributed on the boundary, as demonstrated in Figure 20.For modeling the heat conduction problem with composite materials, the DDM [37] was adopted so that the domain boundary can be divided into several subdomains.Figure 19 shows that the two-dimensional object boundary was divided into several sub-boundaries: Γ 1 , Γ 2 , Γ 3 , . . ., Γ 10 .The Dirichlet boundary condition was imposed on the sub-boundaries using the analytical solution as shown in Equation (52).At the interface, the sub-boundaries, including Γ 2 , Γ 3 , Γ 4 , Γ 5 , Γ 6 , Γ 7 , Γ 8 and Γ 9 , should satisfy the flux conservation and the continuity of temperature.Therefore, the boundary conditions are given as follows: To examine the accuracy of the proposed method, we uniformly collocated 4984 inner points inside the domain to obtain the solution of the temperature.The computed temperature distribution of the inner points is depicted in Figure 21.It was found that the computed results agree very well with those from the analytical solution.To evaluate the accuracy of the proposed method, the absolute error and relative error of the computed results with the analytical solution were calculated, as shown in Figure 22.We obtained highly accurate numerical solutions in the order of 10 −10 and 10 −11 for the absolute error and relative error, respectively.
modeling the heat conduction problem with composite materials, the DDM [37] was adopted so that the domain boundary can be divided into several subdomains.Figure 19 shows that the two-dimensional object boundary was divided into several sub-boundaries: Γ , should satisfy the flux conservation and the continuity of temperature.Therefore, the boundary conditions are given as follows: To examine the accuracy of the proposed method, we uniformly collocated 4984 inner points inside the domain to obtain the solution of the temperature.The computed temperature distribution of the inner points is depicted in Figure 21.It was found that the computed results agree very well with those from the analytical solution.To evaluate the accuracy of the proposed method, the absolute error and relative error of the computed results with the analytical solution were calculated, as shown in Figure 22.We obtained highly accurate numerical solutions in the order of

Discussion
In the CTM, the T-complete functions have to be derived by finding the general solutions of the governing equation.The solutions are expressed as a series using the linear combination of the T-complete functions.Since the T-complete functions satisfy the governing equation, we may place the boundary collocation points only on the boundary of the domain.In addition, the source point of the CTM is only one.The CTM requires the evaluation of a coefficient for each term in the series.On the other hand, the MFS adopts the fundamental solutions as its basis function.Since the fundamental solution has only one term, the MFS has to collocate many source points outside the domain to construct its basis function as a series.The MFS also requires the evaluation of a coefficient for each term in the series.For the CTM, the ill-posedness is often found when higher order terms of the basis functions are used.For the MFS, there is only one basis function but it requires many source points.Although the MFS does not require higher order terms of basis functions, the position of the source points in the MFS, however, is very sensitive to the results due to the singular characteristic of the fundamental solution of the differential operator.
In this article, we propose a novel method using non-singular basis functions instead of the singular solution.The collocation scheme of the proposed method is similar to the MFS.The non-singular basis functions are adopted from the cylindrical harmonics.Due to the adoption of the non-singular basis function, the locations of the source points are not sensitive to the results.This

Discussion
In the CTM, the T-complete functions have to be derived by finding the general solutions of the governing equation.The solutions are expressed as a series using the linear combination of the T-complete functions.Since the T-complete functions satisfy the governing equation, we may place the boundary collocation points only on the boundary of the domain.In addition, the source point of the CTM is only one.The CTM requires the evaluation of a coefficient for each term in the series.On the other hand, the MFS adopts the fundamental solutions as its basis function.Since the fundamental solution has only one term, the MFS has to collocate many source points outside the domain to construct its basis function as a series.The MFS also requires the evaluation of a coefficient for each term in the series.For the CTM, the ill-posedness is often found when higher order terms of the basis functions are used.For the MFS, there is only one basis function but it requires many source points.Although the MFS does not require higher order terms of basis functions, the position of the source points in the MFS, however, is very sensitive to the results due to the singular characteristic of the fundamental solution of the differential operator.
In this article, we propose a novel method using non-singular basis functions instead of the singular solution.The collocation scheme of the proposed method is similar to the MFS.The non-singular basis functions are adopted from the cylindrical harmonics.Due to the adoption of the non-singular basis function, the locations of the source points are not sensitive to the results.This proposed method resolves a major difficulty in the MFS for finding appropriate locations for source points.From the results of the numerical examples, the proposed method was found to be superior in accuracy.
For modeling heterogeneous multi-layer composite materials, the DDM is adopted.The idea of the DDM for modeling composite materials is that the subdomains intersecting at the interface must satisfy flux conservation and the continuity of the temperature for the problem.We successfully integrated the DDM into the hybrid boundary-type meshless method.Two examples of heat conduction problems in heterogeneous multi-layer composite materials were investigated.The results obtained from the examples reveal that the proposed method can be applied to heterogeneous multi-layer composite materials.In addition, the collocation scheme of the proposed hybrid boundary-type meshless method is relatively simple because we only need to place the collocation points on the boundary of the domain.This would be advantageous for engineering problems involving the design of irregular shapes in applied sciences.

Conclusions
This paper presents a study on solving two-dimensional stationary heat conduction problems in composite laminate materials using the novel hybrid boundary-type meshless method.The findings of this study are summarized as follows: 1.
Most applications of the boundary-type meshless method are still limited to homogeneous problems.This study presents pioneering work to investigate the numerical solutions of two-dimensional stationary heat conduction problems in layered heterogeneous composite materials using the novel boundary-type meshless method.We apply the DDM successfully for modeling composite materials in which the domain is split into smaller subdomains.
The subdomains intersecting at the interface must satisfy flux conservation and the continuity of the temperature for the problem.It is found that the proposed method provides a promising solution for solving heat conduction problems with heterogeneity.2.
In the proposed study, we only needed to place collocation points on the boundary.This demonstrated the advantages of the proposed boundary-type meshless method, including the boundary collocation only and high accuracy.Compared to the mesh-based approach, the proposed method is relatively simple and highly accurate.It is therefore advantageous for the analysis of heat conduction problems with a complex shape.Thus, the results obtained show that the proposed method is highly accurate and computationally efficient for modeling heat conduction problems, especially in heterogeneous multi-layer composite materials compared with other meshless methods.Nevertheless, the proposed method is still limited to problems with constant thermal conductivity.Future studies are suggested for problems involving nonlinear behavior.

Figure 1 .
Figure 1.The profile of the composite materials for the analysis.

Figure 2 .
Figure 2. The configuration of the boundary condition.

Figure 3 .
Figure 3. Illustration of collocation points for boundary and source points.

Figure 2 .
Figure 2. The configuration of the boundary condition.

Figure 2 .
Figure 2. The configuration of the boundary condition.

Figure 3 .
Figure 3. Illustration of collocation points for boundary and source points.

Figure 3 .
Figure 3. Illustration of collocation points for boundary and source points.

Figure 4 .
Figure 4. Comparison of computed results of the proposed method.(a) Comparison of results with the analytical solution and (b) comparison of results with the EFG and the IIEFG methods.

Figure 5 .
Figure 5.Comparison of temperature results with the analytical solution along line AB.

Figure 4 .
Figure 4. Comparison of computed results of the proposed method.(a) Comparison of results with the analytical solution and (b) comparison of results with the EFG and the IIEFG methods.

Figure 4 .
Figure 4. Comparison of computed results of the proposed method.(a) Comparison of results with the analytical solution and (b) comparison of results with the EFG and the IIEFG methods.

Figure 5 .
Figure 5.Comparison of temperature results with the analytical solution along line AB.

Figure 5 .
Figure 5.Comparison of temperature results with the analytical solution along line AB.

Figure 6 .
Figure 6.Comparison of the temperature distribution with the analytical solution.

Figure 6 .
Figure 6.Comparison of the temperature distribution with the analytical solution.

24 Figure 7 .
Figure 7.The configuration of the boundary condition.

Figure 8 .
Figure 8. Illustration of the collocation points for boundary and source points.

Figure 7 .
Figure 7.The configuration of the boundary condition.

24 Figure 7 .
Figure 7.The configuration of the boundary condition.

Figure 8 .
Figure 8. Illustration of the collocation points for boundary and source points.Figure 8. Illustration of the collocation points for boundary and source points.

Figure 8 .
Figure 8. Illustration of the collocation points for boundary and source points.Figure 8. Illustration of the collocation points for boundary and source points.

Figure 9 .
Figure 9.Comparison of computed results of the proposed method.(a) Comparison of results with the analytical solution and (b) comparison of results with the EFG and the IIEFG methods.

Figure 10 .
Figure 10.The comparison of the computed results with the analytical solution of temperature along line OA.

Figure 9 .Figure 9 .
Figure 9.Comparison of computed results of the proposed method.(a) Comparison of results with the analytical solution and (b) comparison of results with the EFG and the IIEFG methods.

Figure 10 .
Figure 10.The comparison of the computed results with the analytical solution of temperature along line OA.

Figure 10 .
Figure 10.The comparison of the computed results with the analytical solution of temperature along line OA.

Figure 11 .
Figure 11.Comparison of temperature distribution with the analytical solution.

8 Γ 2 Γ and 6 Γ
. The Dirichlet boundary condition is imposed on the sub-boundaries using the analytical solution as shown in Equation (41) for the problem.At the interface, the sub-boundaries, including , should satisfy the flux conservation and the continuity of temperature.Therefore, the boundary conditions are given as follows:

Figure 11 .
Figure 11.Comparison of temperature distribution with the analytical solution.

Figure 12 .
Figure 12.Illustration of the collocation of boundary and source points.

2 ΩFigure 12 .
Figure 12.Illustration of the collocation of boundary and source points.

Figure 13 .
Figure 13.The comparison of temperature distribution with the analytical solution.

Figure 13 .
Figure 13.The comparison of temperature distribution with the analytical solution.

24 Figure 13 .
Figure 13.The comparison of temperature distribution with the analytical solution.

Figure 14 .
Figure 14.Accuracy of the proposed method (with the comparison of the analytical solution).(a) Absolute error in the first layer, (b) absolute error in the second layer, (c) relative error in the first layer, and (d) relative error in the second layer.

Figure 15 .
Figure 15.The configuration of the boundary condition.

Figure 16 .
Figure 16.Illustration of the collocation of boundary and source points.

5 ΓFigure 15 .
Figure 15.The configuration of the boundary condition.

Figure 15 .
Figure 15.The configuration of the boundary condition.

Figure 16 .
Figure 16.Illustration of the collocation of boundary and source points.

5 ΓFigure 16 .
Figure 16.Illustration of the collocation of boundary and source points.

Figure 17 .Figure 17 .
Figure 17.The comparison of temperature distribution with the analytical solution.

Figure 17 .Figure 18 .Figure 18 .
Figure 17.The comparison of temperature distribution with the analytical solution.

and 2 k
are thermal conductivities in layer 1 and layer 2, respectively.The boundary and source points are uniformly distributed on the boundary, as demonstrated in Figure20.For

Figure 18 .
Figure 18.Comparison of results with the analytical solution.(a) Absolute error in the first layer, (b) absolute error in the second layer, (c) relative error in the first layer, and (d) relative error in the second layer.

10 10 − and 11 10−
for the absolute error and relative error, respectively.

Figure 19 .
Figure 19.The configuration of the boundary condition.Figure 19.The configuration of the boundary condition.

Figure 20 .
Figure 20.Illustration of the collocation of boundary and source points.

Figure 21 .
Figure 21.The comparison of temperature distribution with the analytical solution.

Figure 20 .
Figure 20.Illustration of the collocation of boundary and source points.

24 Figure 20 .
Figure 20.Illustration of the collocation of boundary and source points.

Figure 21 .
Figure 21.The comparison of temperature distribution with the analytical solution.Figure 21.The comparison of temperature distribution with the analytical solution.

Figure 21 .
Figure 21.The comparison of temperature distribution with the analytical solution.Figure 21.The comparison of temperature distribution with the analytical solution.

Figure 22 .
Figure 22.Comparison of results with the analytical solution.(a) Absolute error in the first layer, (b) absolute error in the second layer, (c) relative error in the first layer and (d) relative error in the second layer.

Figure 22 .
Figure 22.Comparison of results with the analytical solution.(a) Absolute error in the first layer, (b) absolute error in the second layer, (c) relative error in the first layer and (d) relative error in the second layer.
Γ D denotes the Dirichlet boundary, and f denotes the boundary where the Neumann boundary condition is given, and Γ N denotes the Neumann boundary.
(17) is a vector of given values from boundary conditions at collocation points with the size of Q × 1. i ≤ Q and j ≤ Q where i is the number of boundary collocation points for Dirichlet boundary condition, j is the number of boundary collocation points for Neumann boundary condition, g 1 , g 2 , ..., g i are the values of the Dirichlet boundary condition, f 1 , f 2 , ..., f j are the values of ∂u/∂n for the Neumann boundary condition.For simplicity, the commercial program MATLAB backslash operator was used to solve Equation(17).

Table 1 .
Comparison of computed results with those from references.

Table 1 .
Comparison of computed results with those from references.

Table 2 .
Comparison of computed results with those from references.

Table 2 .
Comparison of computed results with those from references.