Solutions of Two-Dimensional Heat Transfer Problems by Using Symmetric

C(ξ,x): Symmetric matrice D(ξ,x): A matrice d: Radius of the support domain f: A function FJ: Nodal value of the j th particle G: Normalizing constant g(j): jth particle in the compact support of W(ξ,x) h: The smoothing length K(x) (ξ,x): A matrice L2 : Global error norm M: Total number of particles in the problem domain m: A number N(x): Number of the particles in the compact support domain n: A number P: A matrix for Taylor expansion polynomials P: Source function Q: A matrix for the unknown variables T: Temperature W(ξ,x): Weight function xi,yi,zi: Physical coordinate direction ∂: Derivation operator ∑: Summation symbol λ: The dimensionality of the space ξi: A point in space Δ: The smallest distance between the particle J and its neighboring particles i: Particle index j: Particle index


Introduction
Being the main rival of meshless methods, the finite element method (FEM) is a robust and thoroughly developed method, and it has been widely used in engineering fields due to its versatility for complex geometry and flexibility for many types of linear and nonlinear problems. Most practical engineering problems related to solids and structures are currently solved by using well developed FEM packages that are commercially available. However, the FEM has some inherent shortcomings of numerical methods that rely on the mesh quality and elements. The following limitations of FEM are becoming increasingly evident [1]: 1. High cost in creating an FEM mesh: Creating a mesh for a domain is a prerequisite in using any FEM code. Usually the analyst has to spend most of the time in mesh creation, and it becomes the major component of the cost of a computer aided design (CAD) project especially for three dimensional problems.
2. Low accuracy of stress: Many FEM packages do not predict the stress accurately. The stresses obtained via the FEM are often discontinuous on element boundaries due to the piecewise (or element-wise) continuous nature of the displacement field assumed in the FEM formulation.
3. Difficulty in adaptive analysis: In an adaptive analysis using the FEM, re-meshing (or re-zoning) is required to ensure proper connectivity of elements. To this end, complex, robust and adaptive mesh generation processors have to be developed that are limited to two-dimensional problems. Technical difficulties have precluded the automatic creation of hexahedron meshes for arbitrary three-dimensional domains. In addition, for threedimensional problems, the computational cost of re-meshing at each step is very expensive, even if an adaptive scheme were available. Moreover, an adaptive analysis requires "the mapping" of field variables between meshes in successive stages of the analysis. This mapping process can often lead to additional computation as well as a degradation of accuracy in the solution.  [2] to study three-dimensional (3D) astrophysics problems, has been successfully applied to analyze solid mechanics and transient fluid mechanics problems. However, it has two shortcomings, namely inaccuracy at particles on the boundary and the tensile instability. Many techniques have been developed to alleviate these two deficiencies, among which are the Corrected Smoothed Particle Method (CSPM) [3,4], the Reproducing Kernel Particle Method (RKPM) [5][6][7], the Modified Smoothed Particle Hydrodynamics (MSPH) method [8][9][10][11] and the Symmetric Smoothed Particle Hydrodynamics (SSPH) method [12,13]. The performance of the CSPH, RKPM, MSPH and SSPH in terms of inaccuracy at particles on the boundary and the tensile instability were already discussed in [3][4][5][6][7][8][9][10][11][12]. The MSPH method has been successfully applied to study wave propagation in functionally graded materials, capture the stress field near a crack-tip, and simulate the propagation of multiple cracks in a linear elastic body. The SSPH method is developed to yield symmetric global matrices and has been applied to 2D homogeneous elastic problems successfully.
The SSPH method constructs basis functions that use only locations of particles. These basis functions are found to be similar to those in the Finite Element Methods (FEM) except that the basis for the derivatives of a function need not be obtained by differentiating the function. Needless to say, the basis for the derivatives of a function can be obtained by differentiating the basis function as in the FEM and meshless methods [12].
The SSPH admits a larger class of kernel functions than some other methods such as the SPH, MSPH, RKPM, and the moving least squares (MLS) methods. For finding kernel estimates of derivatives of a function, the SSPH method does not use derivatives of the kernel function while other methods do; instead, the SSPH method uses basis functions different from those employed to approximate the function itself. The kernel function used to generate the basis functions may even be constant.
The SSPH method is applied to some homogeneous solid mechanics problems and proved to be accurate enough to compete with other computational methods. However, originating from its underlying formulations, it has drawbacks for nonhomogeneous problems. Therefore, this study is initiated to reveal the performance of the SSPH method in solving the heat transfer problems, in particular, non-homogeneous problems. It is observed that the SSPH method yields large errors in solving nonhomogeneous problems since it considers nodal values of the forcing term and variation of forcing term in the area among nodes is not considered. Consequently, the error of the SSPH method increases especially if the forcing term is non-smooth. This paper is organized as follows. The SSPH method is described briefly. Then, solutions of different 2D homogeneous and non-homogeneous steady-state heat transfer problems are presented. The conclusions are drawn at the end.

Symmetric Smoothed Particle Hydrodynamics (SSPH)
The value of a function x = (x , x , x ) can be approximated through the finite Taylor Series expansion provided that continuous derivatives up to the order of (n+1) as follow Neglecting the third and higher order terms and introducing the two matrices P(ξ,x) and Q(x), Equation (1) can be written as where T 2 The elements of the matrix Q(x), that are the kernel estimates of a function, its first derivatives and its second derivatives at 1 2 3 x = (x , x , x ) , are known variables to be found from Equation (2). If we multiply both sides of Equation (2) In the compact support of the kernel function W(ξ,x) associated with the point 1 2 3 x = (x , x , x ) , shown in Figure 1, let there be N(x) particles. In the global numbering system, let the particle number of the j th particle in the compact support of W(ξ,x) be g(j). We evaluate Equation (5) at every particle in the compact support of W(ξ,x) and sum both sides of the equation over these particles to arrive at g (2) g (3) g (6) g (7) x g (5) g (4) g (N(x) where ξ g(j) denotes the coordinates of the particle g(j). Equation (6) can be rewritten as where C(ξ,x)=P(ξ,x) T W(ξ,x)P(ξ,x) and D(ξ,x)=P(ξ,x) T W(ξ,x). It is obvious that the matrix C(ξ,x) is symmetric; that is why, this technique is called the SSPH method. The set of simultaneous linear algebraic equations in Equation (7) can be solved for the unknown elements of the matrix Q(x) (Figure 1).
The symmetry of the matrix C(ξ,x) reduces storage requirements and the CPU time needed to solve Equation (7) for Q(x). None of the matrices in Equation (7) involves derivatives of the kernel function. Thus, a much larger class of functions can be used as the kernel function which improves the practicality and usefulness of the method [13].
For the non-singular matrix C(ξ,x), the solution of Equation (7) is given by where K (x) (ξ,x)=C(ξ,x) -1 D(ξ,x). Alternatively, Equation (8) can be written as where F J =f(ξ J ) and M is equal to the total number of particles in the entire domain of interest. The value of the function and its derivatives at the point x are now expressed in terms of the function at all particles in the entire domain [13]. Then, the components of Equation (8) for a 2D problem can be written explicitly as follows

Numerical Examples
In this section, three different boundary value problems of steady state heat transfer in rectangular Cartesian coordinates are solved having the following governing differential equation where h is the smoothing length, λ the dimensionality of the space and G normalizing constant determined by the condition that the integral of the kernel function over the domain is equal to unity. For the circular support domain, the size of the support domain is controlled by the following scaling factor where h is the smoothing length for particle J which is set equal to Δ and Δ is the smallest distance between the particle J and its neighboring particles. 2D homogeneous and non-homogeneous steady-state heat transfer problems are solved for three different uniform particle distributions of 5×10, 9×18 and 17×37. Convergence of the SSPH method for each problem is calculated by using the following global L 2 error norm [14].

Example 1
In this example, we consider that the source function p is zero. Then, the governing equation and essential boundary conditions are given by 17x37 nodes are placed in the problem domain. The smoothing length h is taken as equal to Δ.
When generating the SSPH basis functions, sufficient number of particles should be included in the kernel function's compact support. For the problems requiring the evaluation of 2 nd order derivatives for the collation method, the scaling factor d should be large enough to have at least six particles in the kernel function's compact support. A larger value of d is not recommended because the CPU time required computing the SSPH basis function increases as d increases.
Effects of the scaling factor d on the L2 error norm of temperatures are illustrated in figure 4, where the RGF and RSGF denote Revised Gauss Function and Revised Super Gauss Function, respectively.
It is concluded by examining the curves in figure 4 that the L 2 error norms of temperatures decrease as the number of particles increases and the Revised Super Gauss Function always gives the smallest L 2 error norm of temperatures. When larger values of d are used, the L2 2 error norm of temperatures for different kernel functions is nearly the same.

Example 2
In this example, the source function p is chosen to be non-zero. The governing equation and the essential boundary conditions are given by whose particular solution can be found as follows In solving this example, uniformly distributed 5×10, 9×18 and 17×37 nodes are placed in the problem domain. As in Example 1, the smoothing length h is taken as equal to Δ. The prescribed temperatures on the boundary edges can be found by evaluating Equation (19) on the boundaries. Then, effects of the scaling factor d on the L2 error norms of temperatures are presented in figure 5.
It is clear that the L 2 error norms of temperatures decrease as the number of particles increases. For the 5×10 particle distribution, the L 2 error norms of temperatures for different kernel functions are almost same. The Revised Super Gauss Function always gives the smallest L 2 error norms of temperatures. When the larger value of d is used, the L 2 error norms of temperatures for different kernel functions are nearly the same.

Example 3
In this example, the source function p is selected as an exponential function. The governing equation and the essential boundary conditions are given by Following, uniformly distributed 5×10, 9×18 and 17×37 nodes are placed in the problem domain. As in the previous examples, the smoothing length h is taken as equal to Δ. The prescribed temperatures on the boundary edges can be found by evaluating the analytical solution given by Equation (11). Effects of the scaling factor d on the L 2 error norms of temperatures are illustrated in figure 6.
Showing the convergence of the method, the L 2 error norms of temperatures decrease as the number of particles increases. For the particle distribution of 5×10, the L 2 error norms of temperatures for different kernel functions are almost the same. The Revised Super Gauss Function always gives the smallest L 2 error norms of temperatures. When a larger value of d is used, the L 2 error norms of temperatures for different kernel functions are nearly the same. The convergence of the method by using the Revised Super Gauss Function and Revised Gauss Function is shown in figure 7, where quadratic convergence is observed.

Conclusions
We used the SSPH method to solve 2D homogeneous and nonhomogeneous steady-state heat transfer problems and compared the results obtained by using two different kernel functions and three different particle numbers. The SSPH basis functions are employed to solve three heat transfer problems, two of which are non-homogeneous. According to the numerical results, it is found that the L 2 error norm of temperatures decreases as the number of particles increases, that is the evidence of convergence of the method. However, more particles in the kernel compact support require much more CPU times and also numerical errors increase through the large value of the scaling factor d.
It is observed that optimum value of the scaling factor d to be used for the 2D homogeneous and non-homogeneous steady-state heat transfer problems is 2 and the most suitable kernel function to be used for the SSPH basis function is Revised Super Gauss Function that always yields less error than the Revised Gauss Function. On the other hand, due to the nature of the SSPH formulations that depends on the evaluation of the forcing terms at the nodes, the errors of the SSPH method in nonhomogeneous problems increase regarding to the forcing term, in particular in the existence of non-smooth forcing terms since it does not consider the variation of forcing term in the area among nodes.
Motivated by this fact, future studies will be carried out to develop the performance of the SSPH method in the existence of non-smooth forcing terms, i.e., non-homogeneous problems. Hence, improved accuracy may be obtained for distributed forcing terms.