Strojniški vestnik - Journal of Mechanical Engineering 51(2005)7-8, 476-483
UDK-UDC 536.2
Izvirni znanstveni clanek - Original scientific paper (1.01)
Application of Meshless Methods for Thermal Analysis
Darrell W. Pepper1 and Božidar Šarler2
1Nevada Center for Advanced Computational Methods, University of Nevada Las Vegas, Las Vegas, NV, 89154-402 7, USA, dwpepper@nscee.edu 2Laboratory for Multiphase Processes, Nova Gorica Polytechnic, Nova Gorica, Slovenia
Abstract
Many numerical and analytical schemes exist for solving heat transfer problems. The meshless method is a particularly attractive method that is receiving attention in the engineering and scientific modeling communities. The meshless method is simple, accurate, and requires no polygonalisation. In this study, we focus on the application of meshless methods using radial basis functions (RBFs) – which are simple to implement – for thermal problems. Radial basis functions are the natural generalization of univariate polynomial splines to a multivariate setting that work for arbitrary geometry with high dimensions. RBF functions depend only on the distance from some center point. Using distance functions, RBFs can be easily implemented to model heat transfer in arbitrary dimension or symmetry.
Introduction
For decades, finite difference, finite volume, and finite element methods (FDM/FVM/FEM) have been the dominant numerical schemes employed in most scientific computation. These methods have been used to solve numerous thermal related problems covering a wide range of applications. A common difficulty in these classic numerical methods is the considerable amount of time and effort required to discretize and index domain elements, i.e., creating a mesh. This is often the most time consuming part of the solution process and is far from being fully automated, particularly in 3D. One method for alleviating this difficulty has been to utilize the boundary element method (BEM). The major advantage of the BEM is that only boundary discretization is required rather than domain, thereby reducing the problem by one order. However, the discretization of surfaces in 3-D can still be a complex process even for simple shapes. In addition, these traditional methods are often slowly convergent, frequently requiring the solution of 100’s of thousands of equations in order to get acceptable accuracy.
In recent years, a novel numerical technique called “meshless methods” (or “mesh-free methods”) has been undergoing strong development and has attracted considerable attention from both science and engineering communities. Currently, meshless
methods are now being developed in many research institutions all over the world. Various methods belonging to this family include: Diffuse Element Methods, Smooth Particle Hydrodynamics Methods, Element-Free Galerkin Methods, Partition of Unity Methods, h-p Cloud Methods, Moving Least Squares Methods, Local Petrov-Galerkin Methods, Reproducing Kernel Particle Methods, and Radial Basis Functions.
A common feature of meshless methods is that neither domain nor surface polygonisation is required during the solution process. These methods are designed to handle problems with large deformation, moving boundaries, and complicated geometry. Recently, advances in the development and application of meshless techniques show they can be strong competitors to the more classical FDM/FVM/FEM approaches [1,2], and may likely become a dominant numerical method for solving science and engineering problems in the 21st century. A recent book by Liu [3] discusses meshfree methods, implementation, algorithms, and coding issues for stress-strain problems. Liu [3] also includes Mfree2D, an adaptive stress analysis software package available for free from the web. Atluri and Shen [4] produced a research monograph that describes the meshless method in detail, including much in-depth mathematical basis. They also present comparison results with other schemes.
476
Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483
C         transport variable                                        «           thermal diffusivity (k/pcp)
D         diffusion coefficient                                     e           emissivity
interiorfunctional                                       a           Stefan-Boltzmann constant
boundary functional                                     <P           trial function; field variable
convective film coefficient                           e           dummy variable
thermal conductivity                                    V          approximation function; streamfunction
f
g h
N         number ofnodes"                                    w          vorticity
q         heat flux                                                     Subscripts
Q         heat source/sink                                          i,j         nodal values
r          radial distance                                             I           internal number of node points
t           time                                                           o           initial value; free parameter
T          temperature                                                Superscripts
V        velocity vector                                            n, n+1    known, unknown values
x         horizontal distance                                       -            approximate solution
There exists various types of meshless methods and    conventional  numerical  methods based on mesh
each method has its advantages and disadvantages.    discretizations -thus the ratio of accuracy to CPU is
Intensive research conducted in many major research    likely to be greater for meshless methods
institutions all over the world are now working to       Meshless   methods   hold   promising   alternative
improve the performance of these approaches. In this    approaches for problems involving fluid flow and
study, we focus on the introduction of the basic    heat transfer analyses. The most attractive feature is
concept  of meshless  methods  using  radial basis    the lack of a mesh that is required in the more
functions (RBFs) - whichare simple to implement.         conventional numerical approaches. This becomes
Currently, there are two major approaches in this    particularly interesting in that one can begin to direction: (i) a domain-type meshless method that    conduct adaptive analyses for CFD problems. was developed by Kansa [5] in 1990; (ii) a boundary-type meshless method that has evolved from the    The Meshless Method using RBFs BEM  [6].  Radial basis functions are the natural
generalization of univariate polynomial splines to a    In 1990, Kansa [5] extended the idea of interpolation
multivariate setting. The main advantage of this type    scheme using RBFs to solving various types of
of approximation  is  that  it  works  for  arbitrary    engineering problems. The method is simple and
geometry  with high dimensions and it does not    direct   and   is   becoming   very   popular   in   the
require any mesh. A RBF is a function whose value    engineering community. To illustrate the application
depends only on the distance from some center point.    of the meshless method using Kansa’s method, we
Using   distance   functions,   RBFs   can   be   easily    first consider the elliptic problems. For simplicity,
implemented to reconstruct a plane or surface using    we consider the 2-D Poisson problem with Dirichlet
scattered data in 2-D, 3-D or higher dimensional    boundary condition
spa Similar to  FEM  techniques,   meshless  methods
produce banded system matrices that can be handled     V2T = f(x,y),    (x,y)efi,
in similar fashion. Both sets of methods can utilize     T                                    r                                 (1)
either direct methods based on Gauss elimination or        =g(x,y),        (x,y) gF
CPU time since the creation of shape functions are    [5]  assumed the  approximate  solution could be
more time-consuming and are performed during the    obtained using a linear combination of RBFs computation. However, less time is spent in setting up  computational  nodes.  Results  using  meshless
methods are typically more accurate than                       T(x,y) = ^(rj)                                             (2)
H
Application of meshless methods for thermal analysis
477
Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483
where {T1, T2, …,TN} are the unknown coefficients to be determined, <p(rj) is some form of RBF (trial function), and r is defined as
rj=^(x_xj)2 + (y_yj)2                                (3)
Since multiquadrics (MQ) are infinitely smooth functions, they are often chosen as the trial function for q>, i.e.,
<|)(rj)=>/rj2 + c2=>/(x-xj)2+(y-yj)2+c2
(4)
where c is a shape parameter provided by the user. The optimal value of c is still a subject of research. Other functions such as polyharmonic splines can also be chosen as the trial function.
By direct differentiation of Eq. (6), the first and second derivatives of <p with respect to x and y can be expressed as
d§
9x
9(|>     y - yj
Figure 1. Interior points and boundary points using Kansa’s method.
An implicit time marching scheme can  be used and Eq. (7) becomes
v rj + c2    ^   v rj2 + c2
Tn+1       T n
At
ry2Tn + 1   _  f (T
n
,VTn)
(8)
d2<$>    (y-yj)2 + c2      g>    (x-x)2+c2 (5) dx              . /r2 j_ r-2               dy2
V2                2
r  +c
V2                2
r  +c
Substituting  Eq.  (3)  into  Eq.  (1)  and  using collocation, one obtains
St
(x   -x )   +(y   -y )   +2c
= f(x , y ),
v((xi -xj)2 +(yi -yj)2 +c2)  J
i= 1,2,...,N
where At denotes the time step, superscript n+1 is the unknown (or next time step) value to be solved, and superscript n is the current known value. The approximate solution can be expressed as
T(xi,yi,tn+1)=^Tn+>j(xi,yi)                         (9)
Substituting Eq. (9) into Eq. (8), one obtains
ITj
(x, -x.)  +(y. -y.)  +c   =g(x,y),
i = N +1,N +2,...,N
(6) where NI denotes the total number of interior points and NI+1,…, N are the boundary points. Figure 1 shows two sets of interpolation points: interior and boundary points. Note that Eq. (6) is a linear system of N X N equations and can be solved by direct Gaussian elimination. Once the unknown coefficients {T1, T2, …, TN} are found, the solution of T can be approximated at any point in the domain.
For time dependent problems, we consider the following heat equation as an example:
V^ Tjn+1
At
aV2
(x,y)=— Tn(xi,yi) At
f(xi,yi,tn,Tn(xi,yi),vTn(xi,yi)),    i = 1,2,...,NI
ŁTjn+>(xi,yi)=g(xi,yi,tn+1),    i=NI+1,...,N
(10) which produces an N X N linear system of equations for the unknown Tn+1. Note that the right hand side
of the first equation in Eq. (10) can be updated before the next time step, i.e.,
3T
-aV2T=f(x,y,T,VT)
(7)
478

Pepper D.W. - Šarler B.
Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483
Tn(xi,yi) = XTjn*j(xi,yi),
N
Tn(xi,yi) = XTnj(xi,yi),                           (11)
j=1   j   dx
Tn(xi,yi) = 2JTjn ~(xi , yi) j=1      dy
Heat Transfer Applications
To illustrate the use of meshless methods, let us begin with a simple heat transfer problem. The governing equation for temperature transport can be written as
Table 1. Comparison of results for example 1
Method
cT dt
Exact
V'VT = aV2T+Q
(12)
FEM
BEM
Meshless
mid-pt (oC)   Elements
94.512
94.605
94.471
256
94.514
64
0
Nodes
289
65
325
q + kVT-h(T-T„)-ea(T4-T 4)=0               (13)
T(x,0) = To                                                      (14)
where V is the vector velocity, x is vector space, T(x,t) is temperature, T„ is ambient temperature, To is initial temperature, d is thermal diffusivity (k/Pcp), e is emissivity, ais the Stefan-Boltzmann constant, h is the convective film coefficient, q is heat flux, and Q is heat source/sink. Velocities are assumed to be known and typically obtained from solution of the equations of motion (a separate program is generally used for fluid flow [7]).
In this first example, a two-dimensional plate is subjected to prescribed temperatures applied along each boundary [8], as shown in Fig. 2. The temperature at the mid-point (1,0.5) is used to compare the numerical solutions with the analytical solution. The analytical solution is given as
T-T1       28-,(-1)n1 + 1         nmx     sinh(nry/L)
T-T  nn1
As a second example, a two-dimensional domain is prescribed with Dirichlet and Neumann boundary conditions applied along the boundaries, as shown in Fig. 3(a,b,c). This problem, described in Huang and Usmani [2], was used to assess an h-adaptive FEM technique for accuracy. A fixed temperature of 100oC is set along side AB; a surface convection of 0oC acts along edge BC and DC with h = 750 W/moC and k = 52 W/moC. The temperature at point E is used for comparative purposes. The severe discontinuity in boundary conditions at point B creates a steep temperature gradient between points B and E. Figures 3(b,c) show the initial and final FEM meshes after two adaptations using bilinear triangles. The analytical solution for the temperature at point B is T = 18.2535oC. Table 2 lists the results for the three methods compared with the exact solution. The initial 3-noded triangular mesh began with 25 elements and 19 nodes.
Table 2. Comparison of results for example 2
L     sinh(n7iW / L)
which yields 9(1,0.5) = 0.445, or T(1,0.5) = 94.5oC. Table 1 lists the final temperatures at the mid-point using a finite element method, a boundary element method, and a meshless method, compared with the exact solution.
Method Exact
FEM
BEM
Meshless
Pt E (oC)
18.2535
18.1141
18.2335
18.253 1
Elements	Nodes
0	0
256	155
32	32
0	83
A simple irregular domain is used for the third example and results compared with the three methods. Results from a fine mesh FEM technique (without adaptation) are used as a reference benchmark [7]. The discretized domain and accompanying boundary conditions set along each
Application of meshless methods for thermal analysis

479
Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483
surface are shown in Fig. 4. The FEM results are displayed as contour intervals.
Figure 5(a,b) shows meshless results (using FEM fine mesh nodes for contouring) versus FEM solutions using adapted quadrilateral elements. Heat conduction occurs as a result of constant temperatures set on the top and bottom surfaces, adiabatic faces in the upper right cutout and lower cutout portions, and convective heating along the right and left vertical walls. Adaptive meshing occurs in the corners as a result of steep temperature gradients; this is not evident when using meshless methods.
(c) Figure   3.   Problem   (a)   geometry   -   boundary conditions, (b) initial FEM mesh, and (c) final FEM adapted mesh (from [2])
										
								'    /		
										
										
										
										
				>                                                              L-						1
				q-I						
Figure 4: Problem specification for heat transfer in a user-defined domain.
The FEM, BEM, and meshless mid-point values at (0.5,0.5) are listed in Table 3.
Table 3. Comparison of results for example 3
Method
FEM
BEM
Meshless
mid-pt (oC)
75.899
75.885
75.893
Elements	Nodes 178 37 96
138	
36	
0	
All three techniques provide accurate results for the three example cases. The meshless method was clearly the fastest, simplest, and least storage demanding method to employ. Advances being made in meshless methods will eventually enable the scheme to compete with the FEM and BEM on a much broader range of problems [3,4]. Dr. Y. C. Hon1 is a leading expert in the application of Kansa’s   method.   Much   work   in   engineering
480

Pepper D.W. - Šarler B.
Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483
modeling using Kansa’s method has been done by his    Natural Convection Test Case research group.
Figure 5: FEM solutions (a) meshless (on FEM fine mesh) and (b) adapted mesh.
1 Department of Mechanical Engineering, Hong Kong University, Hong Kong, China
Natural convection within a 2-D rectangular enclosure is a well-known problem commonly used to test the ability of a numerical algorithm to solve for both fluid flow and heat transfer. The equations are strongly coupled through the buoyancy term in the momentum equations and the temperature. There are various ways to nondimensionalize the equations, and numerous references can be found in the literature and on the web regarding these various forms. The solution to the problem generally splits between solving either the primitive equations for velocity or the vorticity equation, coupled with the transport equation for temperature. The issue in this early development of the meshless approach is not to dwell on various schemes to deal with pressure (e.g., projection methods or the SIMPLE scheme both of which are well known). Hence, most researchers that have developed meshless approaches use the streamfunction-vorticity and temperature equations [3]. These equations are the well-known set generally formulated as follows:
 + V-Vco = PrV co- Pr-Ra-VT dt
5T St
+ V-
VT = V2T
V2\|/ = -oo
(20)
(21)
(22)
where a) is vorticity and v is streamfunction, with the conventional definitions for velocity in terms of the streamfunction gradients. Pr is the Prandtl number (v/a) and Ra is the Rayleigh number. Figure 6 shows the physical and computational domain with accompanying boundary conditions. Two types of nodal configurations are shown in Fig. 7 (a,b) utilizing 256 nodes. Results are in excellent agreement with well-known results in the literature for 103< Ra < 105 [3]. Figure 8 (a,b) shows streamlines and isotherms for the differentially heated enclosure for Ra = 105. Convergence rates showing the difference in rates between a conventional FDM and applications of two meshless techniques is discussed in Liu [3]. The two meshless methods converged more rapidly than the finite difference scheme.
Application of meshless methods for thermal analysis

481
Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483
3t
(pC) + V.(pVC)=V.(DVC)+S
(23)
where p, „ V, t, D, and S denote density, transport variable, velocity, time, diffusion matrix, and source. The transport variable C consisted of enthalpy C(h(cp = T)), velocity C(cp = u,v), and pressure C(<p = p), with a pressure correction Poisson equation used to resolve the pressure. The nonlinear equations solved with the meshless technique were of the form
Figure 6. Boundary conditions for natural convection within a rectangular enclosure (from [3]).
v> = e+v.0		(24)
5   , 0=  —(p'C-S) _dt	/D	(25)
0 =[p'VC-D'V(|)l/D		(26)
Figure 7. Nodal configurations for a) uniform distribution and b) arbitrary distribution for 256 nodes (after [3]).
where p' denotes density, C(ch ) the transport variable, t is time, V is velocity, and D is te diffusion matrix with D' being the nonlinear anisotropic part. The variable C(<p) denotes the relation between the transported and the diffused variable. The solution requires the use of an iterative technique. The final form of the transformed Poisson equation is
V2(|)= 0 + 0,,,
(c|>-(|>) + V'0 + V"04((|>-(|>)
(27)
where the bar denotes values from the previous iteration. Time discretization utilizes the relation
"pC(^)-pC(^o) 0«  --------------------S
At
/ D
(28)
Figure 8. Natural convection results showing a) streamlines and b) isotherms for Ra = 105 using the MLPG method (after [3]).
Sarler et al [9] simulated natural convection within a rectangular enclosure using the RBF approach of Kansa [5]. Solving a nonlinear Poisson re-formulation of the general transport equation representing mass, energy, and momentum, the problem was solved by dividing the physical domain into two parts consisting of an internal array of nodes and a set of boundary nodes for the Dirichlet and Neumann conditions. The governing equation for the transport variable is of the form (with C(q>))
with the unknown field <p approximated by the N global approximation functions Ąn(p) and their coefficients ?n, i.e.
^(p)-Vn(p)qn,    n = 1,2,...,Nr                             (29)
The global radial basis function approximation was based on multiquadrics with the free parameter ro:
Vn

(30)
The coefficients were calculated from the N collocation equations of which Nr were equally distributed over boundary T and Nn over the domain Q. Separate relations were established for the boundary condition indicators.
482

Pepper D.W. - Šarler B.
Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483
The computational domain was discretized into 80 boundary nodes and 361 domain nodes. The multiquadrics constant ro was set to 0.2. Steady state results were achieved after 34 iterations for Ra = 103, 187 iterations for Ra = 104, and 293 iterations for Ra = 105. The calculated values for temperature and velocity were in excellent agreement with results obtained using a fine grid FDM [10].
Conclusions
Meshless methods are a unique and novel numerical technique now making inroads into various fields. Their advantages in solving problems associated with crack propagation and stress/strain including deformation over more conventional numerical schemes have been demonstrated repeatedly in the literature. The application of meshless methods for heat transfer is equally advantageous; such methods have become very competitive with both finite volume and finite element methods for problems involving irregular geometries. The requirements for creating grids as well as the detailed input necessary for establishing volume or element properties is greatly reduced or eliminated. However, much has yet to be done before meshless methods can handle a wide range of fluid flow problems and produce results with confidence and surety. Only a small portion of incompressible flow problems has been addressed; advances are just now being made in the area of porous media flows. Application to compressible flows has yet to be addressed.
While meshless methods may be more accurate than FDM/FVM/FEM techniques, they can be much slower with regards to computational time to achieve convergence. This is due in part to some of the effort needed for numerical integration and subsequent use of a direct matrix solver. However, meshless methods do not need any prior knowledge of their nodal arrangement, as in conventional numerical schemes. This makes the method particularly attractive for developing adaptive capabilities. Since much of a modeler’s efforts are generally spent on developing a good mesh that will lead to a converged solution, the overall time for obtaining problem solutions using meshless methods can be significantly less.
Additional information regarding FDM, FVM, and FEM algorithms and some of the meshless techniques can be obtained from the web site http://www.unlv.edu/NCACM.
Acknowledgements
We wish to thank Professor C. S. Chen and Professor Jichun Li from the Department of Mathematics at UNLV for their helpful insight and assistance.
References
1. R. W. Lewis, K. Morgan, H. R. Thomas, and K. N. Seetharamu, The Finite Element Method in Heat Transfer Analysis, J. Wiley & Sons: Chichester, UK (1996).
2.H-C. Huang and A. S. Usmani, Finite Element Analysis for Heat Transfer, Springer-Verlag: London, UK (1994).
3. G. R. Liu, Mesh Free Methods: Moving Beyond the Finite Element Method, CRC Press, Boca Raton, FL (2002).
4. S. N. Atluri and S. Shen, The Meshless Local Petrov Galerkin (MLPG) Method, Tech Science Press, Encino, CA (2002).
5.E. J. Kansa, Multiquatric – A Scattered Data Approximation Scheme with Applications to Computational Fluid Dynamics II, Computers Math. Appl., 19, 8/9, (1990), 147-161.
6. S. N. Atluri, H. K. Kim, and J. Y. Cho, J.Y., A Critical Assessment of the Truly Meshless Local Petrov-Galerkin (MLPG), and Local Boundary Integral Equation (LBIE) Methods, Computational Mechanics, 24, (1999), 348-372.
7. D. W. Pepper, D. B. Carrington, and L. Gewali, A Web-based, Adaptive Finite Element Scheme for Heat Transfer and Fluid Flow. ISHMT/ASME 4th Conf. on Heat and Mass Transfer, Jan. 12-14, Pune, India (2000).
8. F. P. Incropera and D. P. DeWitt, Fundamentals of Heat and Mass Transfer, 5th Ed., J. Wiley & Sons: New York (2002).
9.B. Sarler, J. Perko, C. S. Chen, and G.Kuhn, A Meshless Approach to Natural Convection, S. N. Atluri and D. W. Pepper (Eds.), Proceedings of the 22nd ICES Conference, July 31-Aug. 2, 2002, Reno, NV (2002).
10.D. W. Pepper and C. S. Chen, A Meshless Method for Modeling Heat Transfer, S. N. Atluri and D. W. Pepper (Eds.), Proceedings of the 22nd ICES Conference, July 31-Aug. 2, 2002, Reno, NV (2002).
Application of meshless methods for thermal analysis                                           483