A cell-based smoothed radial point interpolation method (CS-RPIM) for heat transfer analysis

Engineering Analysis with Boundary Elements 40 (2014) 147–153

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

A cell-based smoothed radial point interpolation method (CS-RPIM) for heat transfer analysis X.Y. Cui a,b, S.Z. Feng a, G.Y. Li a,n a b

State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha 410082, PR China State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, PR China

art ic l e i nf o

a b s t r a c t

Article history: Received 13 July 2013 Accepted 3 December 2013 Available online 9 January 2014

A cell-based smoothed radial point interpolation method (CS-RPIM) is further extended to solve 2D and 3D heat transfer problems. For this method, the problem domain is first discretized using triangular elements or tetrahedral elements, and each element is further divided into several smoothing cells. Then, the field functions are approximated using RPIM shape functions. Finally, the CS-RPIM utilizes the smoothed Galerkin weak form to obtain the discretized system equations in these smoothing cells. Several numerical examples with different kinds of boundary conditions are investigated to verify the validity of the present method. It has been found that the CS-RPIM can achieve better accuracy and higher convergence rate, when dealing with the 2D and 3D heat transfer analysis. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Heat transfer RPIM CS-RPIM Meshfree method

1. Introduction Heat transfer analysis [1–4] is of great importance in many practical engineering problems. As it is quite difficult to find analytical solution for such problems with complex geometry and boundary conditions, numerical methods are widely employed to deal with the related analysis. In the past several decades, various meshfree methods have been developed and extended to deal with these problems. For instance, Liu [5] employed the meshless local Petrov– Galerkin (MLPG) method to deal with the heat transfer analysis. Singh [6] utilized the element free Galerkin (EFG) method to obtain the numerical solution of heat transfer problems, in which the approximation function is constructed from a set of scattered nodes. Khosravifard [7] proposed an improved meshless radial point interpolation method (PRIM) for the analysis of nonlinear transient heat conduction problems. Recently, a node-based smoothed point interpolation method (NS-PIM) is applied by Wu [8,9] to analyze heat transfer process, which can produce an upper bound solution. Wu [10] also extend an edge-based smoothed point interpolation method (ES-PIM) to analyze 2D and 3D transient heat transfer problems, which can provide more accurate solutions than standard FEM when using the same mesh. Incorporating meshfree techniques with the standard FEM, Liu [11–13] proposed a smoothed finite element method (SFEM) by using the strain smoothing technique in FEM settings. The SFEM

n

Corresponding author. Tel.: þ 86 731 8821717; fax: þ 86 731 8822051. E-mail address: [email protected] (G.Y. Li).

0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.12.004

further divides elements into several smoothing cells and computes the integrals along the edge of these smoothing cells. Using the SFEM idea, a cell-based smoothed radial point interpolation method (CS-RPIM) based on the generalized gradient smoothing operation has been proposed for static and free vibration analysis of 2D solids [14]. In this work, 2D and 3D heat transfer problems are solved by the CS-RPIM. For this method, the computational domain is first discretized into a set of triangular cells or tetrahedral cells, and each cell is further divided into smoothing cells. The discretized system equations are derived using the smoothed Galerkin weak form. Numerical examples with various kinds of boundary conditions are presented to illustrate the validity of the CS-RPIM for the heat transfer analysis.

2. Thermal governing equations and boundary conditions The differential equations governing heat conduction are given as ! ∂Tðx; y; z; tÞ ∂2 T ∂2 T ∂2 T ρc ¼ kx 2 þ ky 2 þ kz 2 þQ ðx; y; z; tÞ ð1Þ ∂t ∂x ∂y ∂z where Tðx; y; z; tÞ is the temperature at time t, Q ðx; y; z; tÞ is the rate of internal heat generation, ρ is the density, c is the specific heat, kx, ky and kz are the thermal conductivities in the x, y, z directions, respectively. The initial condition and thermal boundary conditions are simply stated here as Initial condition : T ¼ T 0

ð2Þ

148

X.Y. Cui et al. / Engineering Analysis with Boundary Elements 40 (2014) 147–153

Dirichlet boundary : T j Γ ¼ T w

ð3Þ

 ∂T  ¼q Neumann boundary :  k ∂n0 Γ

ð4Þ

n

 ∂T  Robin boundary :  k ¼ hc ðT  T 1 Þ ∂n0 Γ Adiabatic boundary :  k

A constraint condition are given here [16]

 ∂T  ¼0 ∂nΓ

ð5Þ

ð6Þ

where T 0 is the initial temperature, T w is the known temperature, T 1 is the environmental temperature, q is the prescribed heat flux, hc is the convection coefficient, n0 is the unit outward normal to the boundary and Γrepresents the boundary.

∑ pj ðxk Þai ¼ 0;

i¼1

Combining Eqs. (10) and (14) yields the following equations: #   "R   Pm  a  q Ts a ~ ¼ ¼G Ts ¼ ð15Þ T Pm 0 0 b b Solving Eq. (15) yields        Us Ts a a ¼ G1 ¼ G1 0 0 b b The approximation of function T is finally expressed as   h i Ts T T ¼ φTs T ¼ Rq Pm G  1 0

φ ¼ ½ φ1 This section presents the RPIM approximation for the temperature field using local nodes and the radial basis functions (RBF) augmented with polynomial basis functions [15]. Consider a problem domain with a set of arbitrarily scattered points xi, (i¼ 1, 2,…, n ), n is the number of nodes in the local support domain. The approximation of temperature function T can be expressed as m

i¼1

j¼1

φ2

⋯

φn 

ð7Þ

where Ri(x) is radial basis function, ai is the unknown coefficient for functions Ri(x), bj is the coefficient for polynomial basis Pj(x), n is the number of field nodes in the local support domain and m is determined according to the polynomial basis. In this work, the multi-quadrics RBF (MQ-RBF) is utilized, which is given as Ri ðxÞ ¼ ½ðx  xi Þ2 þ ðy yi Þ2 þ ðz  zi Þ2 þ ðαc dc Þ2 q

ð8Þ

where dc is the equivalent length of the background cell, q and αc are shape parameters, which are real and arbitrary. This characteristic has been examined in detail by Liu [15]. The polynomial basis function for two dimensional domains is given as PT ðxÞ ¼ ½1; x; y; z; …

ð9Þ

The coefficients in Eq. (7) can be determined by enforcing the field function to be satisfied at the n nodes within the local support domain of the point of interest x. This measure leads to n linear equations, which can be expressed in the matrix form as Ts ¼ R q a þ Pm b

ð16Þ

ð17Þ

ð18Þ

where n

m

i¼1

j¼1

φk ðxÞ ¼ ∑ Ri ðxÞGði;kÞ þ ∑ pj ðxÞGðn þ j;kÞ

ð19Þ

where Gði;kÞ is the element of matrix G  1 . So, the approximation function T can be given as n

T ¼ ∑ Ri ðxÞai þ ∑ pj ðxÞbj ¼ RT ðxÞa þ PT ðxÞb

ð14Þ

where φðxÞ is RPIM shape functions, which is given by

3. Radial point interpolation method

n

j ¼ 1; 2; …; m

T ¼ ∑ φi T i

ð20Þ

i¼1

4. Heat transfer analysis using CS-RPIM 4.1. 2D heat transfer analysis In the 2D analysis, the problem domain is first discretized using triangular elements as in the FEM. The triangular cell is termed as “parent” cell for convenience in our discussion. The edge of a parent cell is called “cell edge”. Each “parent” cell is further divided into non-overlapping smoothing cells Ωk ¼ [ SC C ¼ 1 ΩkðCÞ , as shown in Fig. 1. The edge of the smoothing cell is called “segment”. All the segments of the smoothing cell form the boundary of the smoothing cell, and the boundary is denoted as Γk(C). The temperature gradient in the Cth smoothing cell Ωk(C) in Ωk is assumed constant which is expressed using following generalized gradient smoothing

ð10Þ

where Ts is the vector of function values, Ts ¼ fT 1 ; T 2 ; …; T n gT

ð11Þ

Rq is the moment matrix of RBFs, which can be given as 2 3 R1 ðx1 Þ R2 ðx1 Þ ⋯ Rn ðx1 Þ 6 7 6 R1 ðx2 Þ R2 ðx2 Þ ⋯ Rn ðx2 Þ 7 7 Rq ¼ 6 6⋯ 7 ⋯ ⋯ ⋯ 4 5 R1 ðxn Þ R2 ðxn Þ ⋯ Rn ðxn Þ

ð12Þ

nn

and matrix Pm 2 p1 ðx1 Þ 6 6 p1 ðx2 Þ Pm ¼ 6 6⋮ 4 p1 ðxn Þ

is defined as p2 ðx1 Þ

⋯

p2 ðx2 Þ

⋯

⋮

⋮

p2 ðxn Þ

⋯

pm ðx1 Þ

3

7 pm ðx2 Þ 7 7 7 ⋮ 5 pm ðxn Þ

ð13Þ

Fig. 1. The problem domain is divided into triangular background cells for 2D analysis. Each triangular cell is called parent cell. The kth parent cell is further divided into SC smoothing cells.

X.Y. Cui et al. / Engineering Analysis with Boundary Elements 40 (2014) 147–153

149

3 Cell I

Field node Centroid of the surface triangle Centroid of the tetrahedron

Node k

Mid-edge-point

Fig. 3. Illustration of constructing smoothing domains for 3D problems.

1

2 SC=3

Fig. 2. A parent cell is divided into 3 smoothing cells. The circle presents the field node, and the square denotes the integration sampling point.

technique [17]: Z 1 εTkðCÞ  n UT dΓ AkðCÞ ΓkðCÞ

ð21Þ

where T is the temperature field in the smoothing cell Ωk(C), Ak(C) is the area of the Cth smoothing cell, and n is the outward normal matrix containing the components of the outward normal vector on the boundary Γk(C). It has been demonstrated that the SC¼3 scheme always gives more accurate results in the 2D static and free vibration analysis [14], so triangular background cell with 3 smoothing cells is considered in this work, as shown in Fig. 2. Substituting Eq. (20) into Eq. (21), the smoothed temperature gradient can be rewritten as ! Z NP 1 NP εTkðCÞ  n Uφi dΓ Ti ¼ ∑ gikðCÞ Ti ð22Þ ∑ AkðCÞ i ¼ 1 Γ kðCÞ i¼1 where NP is the number of the supporting nodes of point x, Ti is the temperature at node i, and gikðCÞ is the smoothed temperature gradient interpolation matrix, which is given by 2 3 i 6 bkxðCÞ 7 i ð23Þ gkðCÞ ¼ 4 i 5 bkyðCÞ in which i bkxðCÞ i bkyðCÞ

¼ ¼

1

N sgem

AkðCÞ

I¼1

1

N sgem

AkðCÞ

I¼1

∑

∑

"

NG

#

nIx lI ∑ φi ðxIJ Þ "

J¼1 NG

nIy lI ∑ φi ðxIJ Þ

# ð24Þ

J¼1

where lI is the length of Ith segment of the smoothing cell, Nsgem is the number of the segment of the smoothing cell, NG is the number of the Gauss points used in each segment, nIx and nIy are the components of the outward unit normal to the Ith boundary segment, and xIJ is the coordinate vector of the Jth gauss point on the Ith segment.

cells, as shown in Fig. 3. The face of the smoothing cell is called “segment”. All the segments of the smoothing cell form the boundary of the smoothing cell, and the boundary is denoted as Ωk(C). The temperature gradient in the Cth smoothing cell is also calculated using the generalized gradient smoothing technique, which is similar to the 2D analysis: ! Z NP 1 NP T εkðCÞ  n Uφi dΩ Ti ¼ ∑ gikðCÞ Ti ð25Þ ∑ V kðCÞ i ¼ 1 ΩkðCÞ i¼1 where Vk(C) is the volume of the Cth smoothing cell, n is the outward normal matrix on the boundary Ωk(C), NP is the number of the supporting nodes of point x, Ti is the temperature at node i, and gikðCÞ is the smoothed temperature gradient interpolation matrix for the 3D heat transfer problems, which is given by 2 i 3 bkxðCÞ 6 7 6 i 7 7 ð26Þ gikðCÞ ¼ 6 6 bkyðCÞ 7 4 i 5 bkzðCÞ in which

i bkyðCÞ i bkzðCÞ

In the 3D analysis, the problem domain is discretized using tetrahedral elements. The tetrahedral cell is also termed as “parent” cell. The face of a parent cell is called “cell face”. Each “parent” cell is further divided into non-overlapping smoothing

¼

"

∑

V kðCÞ

I¼1

1

N sgem

V kðCÞ

I¼1

1

Nsgem

V kðCÞ

I¼1

NG

#

nIx AI ∑ φi ðXIJ Þ J¼1

"

NG

#

nIy AI ∑ φi ðXIJ Þ

∑

J¼1

"

NG

nIz AI ∑ φi ðXIJ Þ

∑

# ð27Þ

J¼1

Considering the accuracy and stability, tetrahedral background cell with 4 smoothing cells is utilized in this work. After obtaining the smoothed temperature gradient, the smoothed Galerkin weak form is utilized to deal with Eq. (1). Then, the discretized system equations can be obtained     n o ∂T ð28Þ ¼ Qc þ Qq KT fTg þ NT ∂t where

2

Z KT ¼

Ω

Z NT ¼

4.2. 3D heat transfer analysis

¼

N sgem

1

i

bkxðCÞ ¼

Ω

kx

6 gT 4 0 0

0 ky 0

0

3

0 7 5g dV kz

ρcHT H dV

ð29Þ

ð30Þ

Z fQ c g ¼

Γ

HT hc ðT 1  TÞdΩ

ð31Þ

qHT dΩ

ð32Þ

Z fQ q g ¼

Γ

150

X.Y. Cui et al. / Engineering Analysis with Boundary Elements 40 (2014) 147–153

where KT is the smoothed conduction matrix, NT is the thermal capacitance matrix, fQ c g, fQ q g are the thermal load vectors caused by different boundary conditions which have been represented previously. The backward difference technique is used for time approximation, which can avoid the oscillation especially with longer time step. Therefore, f∂TðtÞ=∂tg is expressed as   ∂TðtÞ fTðtÞg  fTðt  ΔtÞg ¼ þ ΟðΔtÞ ð33Þ ∂t Δt Substituting Eq. (33) into Eq. (28), the following equation can be obtained.        n o  NT  NT  TðtÞ ¼ Q c þ Q q þ Tðt  ΔtÞ ð34Þ KT þ Δt Δt Once the initial temperature field is given, the temperature field at arbitrary time t can be obtained by Eq. (34).

5. Results and discussion Three numerical examples are investigated to validate the stability, accuracy, and convergence rate of CS-RPIM, when dealing with 2D and 3D heat transfer problems. For comparison, the FEM codes are also developed to evaluate the same problems using the same meshes as the CS-RPIM. To study the convergence, the thermal equivalent energy is defined as 2 3 kx 0 0 Z 6 7 ð35Þ U T ¼ gT 4 0 ky 0 5g dΩ Ω 0 0 kz As the analytical solutions of these complex problems are not available, reference solution is obtained using the finite element commercial software ABAQUS, in which very fine meshes with high-order elements are utilized. It should be noted that the RPIM shape functions are created using the MQ-RBF augmented with the linear polynomial basis (m ¼3), and the shape parameters q is taken as 1.03, αc is taken as 0.35. The reason for such choices can be found in [18].

hc ¼100 W/m2 1C, T 1 ¼200 1C. As the thermal conductivity is considered to vary with temperature, the Newton–Raphson iteration method is utilized to deal with the nonlinearity. The CS-RPIM and FEM solutions are obtained using the same triangular meshes with 231 nodes. The reference solution is obtained using ABAQUS quadrilateral elements with a large number of nodes (11 896) for comparison purposes. In addition, coordinates of the irregular mesh are also generated to study the sensitivity of CS-RPIM for different meshes using the following equations: xir ¼ x þ Δx Ur c U αir yir ¼ y þΔy U r c Uαir

where Δx and Δy are the initial regular nodal spacing in x and y direction, respectively. rc is a computer generated random number between  1.0 and 1.0, αir is a prescribed irregularity factor chosen between 0.0 and 0.5. Irregular meshes with αir ¼ 0.1 and αir ¼0.5 are shown in Fig. 5. Fig. 6 shows the comparison of temperature distributions along the central line CD, which are obtained by CS-RPIM, FEM, and ABAQUS, respectively. It can be clearly seen that the results of CSRPIM agree well with the reference solutions, and present higher precision than the FEM when using the same triangular meshes. This kind of improved accuracy is mainly due to the soften effect caused by the smoothing operation. To study the stability of present method, Figs. 7 and 8 further show the nodal temperature of point A and B obtained by different methods using irregularly distributed field nodes. It is clearly shown that the present CSRPIM are insensitive to the mesh distortion even when the meshes are severely distorts with αir ¼0.5.

5.1. A rectangular plate To verify the present CS-RPIM formulation, 2D nonlinear steady heat transfer problem is first studied as illustrated in Fig. 4. In the computation, the related parameters are taken as kx, ky ¼ (10þ0.1T) W/m 1C, l1 ¼1.2 m, l2 ¼ 2.0 m. For the Neumann boundary, q¼ 1000 W/m2, and for the Robin boundary

y l2 B

Neumann boundary C

D l1

Robin boundary

A

Adiabatic

o

x Adiabatic Fig. 4. Geometry and the boundary conditions of the rectangular plate.

ð36Þ

Fig. 5. Irregular meshes with different irregularity factors αir.

X.Y. Cui et al. / Engineering Analysis with Boundary Elements 40 (2014) 147–153

Fig. 6. Comparison of temperature distribution along the CD line (C-D).

151

Fig. 8. Temperature of point B for different methods using irregularly distributed field nodes.

A

Robin boundary B

Fig. 7. Temperature of point A for different methods using irregularly distributed field nodes.

5.2. A square plate with a cylindrical hole This section represents 3D transient heat transfer analysis of a square plate with the dimensions in x–y plane as [  0.1,0.1]  [  0.1,0.1] and the thickness in z-direction as [0, 0.02]. The plate has a central cylindrical hole in z-direction with radius r ¼ 0.05 m. Due to the symmetric property, only a quarter of the plate is modeled. In addition, the plate is subjected to the Robin and Dirichlet boundary conditions, as shown in Fig. 9. Some computational parameters are given as kx, ky, kz ¼15.0 W/m 1C, T0 ¼ 50 1C, ρ¼7800 kg/m3 and c ¼125 J/kg 1C. For the Robin boundary on the upper plane, T 1 ¼100 1C and hc ¼200 W/m2 1C and for the Dirichlet boundary, T w ¼200 1C. The problem domain is discretized into tetrahedral meshes with 354 nodes, and a reference solution is also obtained using ABAQUS hexahedral elements with a number of nodes (18 312) for comparison. Figs. 10 and 11 plot the nodal temperature history of points A and B calculated by CS-RPIM, FEM, and ABAQUS, respectively. Fig. 12 shows the convergence of thermal equivalent energy for a range of meshes, when the system has already reached its steady state. It can be clearly observed that the results obtained by CSRPIM agree better with reference solutions than those obtained by

r Dirichlet z x boundary y Fig. 9. Schematic of a square plate with mixed boundary conditions.

FEM and the CS-RPIM can achieve much higher convergence rate when dealing with the 3D transient heat transfer analyses. 5.3. An engine pedestal In this section, an engine pedestal is investigated, which is subjected to the Neumann and Robin boundary conditions. The geometry and the boundary conditions are shown in Fig. 13. The related parameters used in the computation are taken as kx, ky, kz ¼25.0 W/m 1C, l1 ¼ l2 ¼0.14 m, r1 ¼0.06 m, r2 ¼0.04 m, R¼ 0.03 m, h1 ¼0.015 m and h2 ¼0.05 m, ρ¼1600 kg/m3 and c¼200 J/kg 1C. For the Neumann boundary, q¼  6000 W/m2, and for the Robin boundary hc ¼300 W/m2 1C, T 1 ¼ 200 1C. This problem is analyzed using the CS-RPIM and FEM with the same tetrahedral mesh (606 nodes). In addition, the reference solution is obtained using ABAQUS with very dense quadrilateral mesh (20 583 nodes) for comparison.

152

X.Y. Cui et al. / Engineering Analysis with Boundary Elements 40 (2014) 147–153

Fig. 13. Schematic of the engine pedestal with Robin and Neumann boundary conditions. Fig. 10. Comparison of temperature history at reference point A.

FEM (606 nodes)

CS-RPIM (606 nodes)

Fig. 11. Comparison of temperature history at reference point B.

Reference (20583 nodes)

Fig. 14. Comparison of the temperature fields.

Fig. 12. Convergence of the thermal equivalent energy versus the number of degrees of freedom.

CS-RPIM perform much better than the standard FEM when dealing with the 3D transient heat transfer problems.

6. Conclusion Fig. 14 shows the comparison of temperature field, when the system has already reached its steady state. Figs. 15–17 further show the comparison of temperature gradient in the x, y, z direction, respectively. These results again verify that the

In this paper, a cell-based smoothing radial point interpolation method (CS-RPIM) is further extended to solve 2D and 3D heat transfer problems. Formulations based on the smoothed Galerkin

X.Y. Cui et al. / Engineering Analysis with Boundary Elements 40 (2014) 147–153

3. When dealing with the heat transfer analysis, the CS-RPIM can achieve much higher convergence rate than the standard FEM using the same meshes.

1200 1000

Temperature gradient ( /m)

153

800 600 400 200 0

Acknowledgments

-200 -400 -600 -800 -1000 -1200

Reference

CS-RPIM

FEM

Fig. 15. Comparison of the temperature gradient in the x direction. 1400

The support of National Science Foundation of China (11002053), Project supported by the State Key Program of National Natural Science of China (61232014), National Basic Research Program of China (2010CB328005), State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology (GZ1212) and Science Fund of State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body (51075005) are gratefully acknowledged.

1200

Temperature gradient ( /m)

1000

References

800 600 400 200 0 -200 -400 -600 -800 -1000 -1200 -1400

Reference

CS-RPIM

FEM

Fig. 16. Comparison of the temperature gradient in the y direction. 600

Temperature gradient ( /m)

400 200 0 -200 -400 -600 -800 -1000 -1200 -1400 -1600

Reference

CS-RPIM

FEM

Fig. 17. Comparison of the temperature gradient in the z direction.

weak form are first presented, and numerical examples with different kinds of boundary conditions are then investigated to examine the stability, accuracy, and convergence rate of the present method. From these studies, several conclusions can be drawn as follows: 1. In the CS-RPIM, there is no need to obtain the derivatives of meshfree shape functions. 2. The CS-RPIM models are very stable and accurate, even when the nodes are extremely distributed.

[1] Islamoglu Y, Kurt A. Heat transfer analysis using ANNs with experimental data for air flowing in corrugated channels. Int J Heat Mass Transf 2004;47(6–7):1361–5. [2] Kim MY. A heat transfer model for the analysis of transient heating of the slab in a direct-fired walking beam type reheating furnace. Int J Heat Mass Transf 2007;50(19–20):3740–8. [3] Naphon P, Kornkumjayrit K. Numerical analysis on the fluid flow and heat transfer in the channel with V-shaped wavy lower plate. Int Commun Heat Mass Transf 2008;35(7):839–43. [4] Aprin L, Mercier P, Tadrist L. Local heat transfer analysis for boiling of hydrocarbons in complex geometries: a new approach for heat transfer prediction in staggered tube bundle. Int J Heat Mass Transf 2011;54(19–20): 4203–19. [5] Liu LH, Tan JY. Meshless local Petrov–Galerkin approach for coupled radiative and conductive heat transfer. Int J Therm Sci 2007;46(7):672–81. [6] Singh A, Singh IV, Prakash R. Meshless element free Galerkin method for unsteady nonlinear heat transfer problems. Int J Heat Mass Transf 2007;50(5–6):1212–9. [7] Khosravifard A, Hematiyan MR, Marin L. Nonlinear transient heat conduction analysis of functionally graded materials in the presence of heat sources using an improved meshless radial point interpolation method. Appl Math Modell 2011;35(9):4157–74. [8] Wu SC, Liu GR, Zhang HO, Zhang GY. A node-based smoothed point interpolation method (NS-PIM) for thermoelastic problems with solution bounds. Int J Heat Mass Transf 2009;52(5–6):1464–71. [9] Wu SC, Liu GR, Zhang HO, Xu X, Li ZR. A node-based smoothed point interpolation method (NS-PIM) for three-dimensional heat transfer problems. Int J Therm Sci 2009;48(7):1367–76. [10] Wu S, Liu G, Cui X, Nguyen T, Zhang G. An edge-based smoothed point interpolation method (ES-PIM) for heat transfer analysis of rapid manufacturing system. Int J Heat Mass Transf 2010;53(9):1938–50. [11] Dai KY, Liu GR, Nguyen TT. An n-sided polygonal smoothed finite element method (nSFEM) for solid mechanics. Finite Elem Anal Des 2007;43(11–12): 847–60. [12] Dai KY, Liu GR. Free and forced vibration analysis using the smoothed finite element method (SFEM). J Sound Vib 2007;301(3–5):803–20. [13] Liu G, Nguyen T, Dai K, Lam K. Theoretical aspects of the smoothed finite element method (SFEM). Int J Numer Methods Eng 2007;71(8):902–30. [14] Cui X, Liu G, Li G. A cell-based smoothed radial point interpolation method (CS-RPIM) for static and free vibration of solids. Eng Anal Bound Elem 2010;34 (2):144–57. [15] Liu G-R. Meshfree methods: moving beyond the finite element method. CRC Press Llc; 2010. [16] Golberg MA, Chen CS, Bowman H. Some recent results and proposals for the use of radial basis functions in the BEM. Eng Anal Bound Elem 1999;23 (4):285–96. [17] Liu G. A generalized gradient smoothing technique and the smoothed bilinear form for Galerkin formulation of a wide class of computational methods. Int J Comput Methods 2008;5(02):199–236. [18] Wang J, Liu G. On the optimal shape parameters of radial basis functions used for 2-D meshless methods. Comput Meth Appl Mech Eng 2002;191(23): 2611–30.