† Corresponding author. E-mail:

Kernel gradient free-smoothed particle hydrodynamics (KGF-SPH) is a modified smoothed particle hydrodynamics (SPH) method which has higher precision than the conventional SPH. However, the Laplacian in KGF-SPH is approximated by the two-pass model which increases computational cost. A new kind of discretization scheme for the Laplacian is proposed in this paper, then a method with higher precision and better stability, called Improved KGF-SPH, is developed by modifying KGF-SPH with this new Laplacian model. One-dimensional (1D) and two-dimensional (2D) heat conduction problems are used to test the precision and stability of the Improved KGF-SPH. The numerical results demonstrate that the Improved KGF-SPH is more accurate than SPH, and stabler than KGF-SPH. Natural convection in a closed square cavity at different Rayleigh numbers are modeled by the Improved KGF-SPH with shifting particle position, and the Improved KGF-SPH results are presented in comparison with those of SPH and finite volume method (FVM). The numerical results demonstrate that the Improved KGF-SPH is a more accurate method to study and model the heat transfer problems.

Grid-based numerical methods such as the finite element method (FEM) and finite difference method (FDM) have been widely applied to various fields of computational fluid dynamics (CFD) and computational solid dynamics (CSD). However, conventional grid-based methods suffer some inherent difficulties in many aspects, such as mesh generation for irregular or complex geometry and cases with heavily distorted mesh. So many meshless methods have been developed over the past twenty years, such as element-free Galerkin method,^{[1–3]} complex variable meshless method,^{[4,5]} finite point method,^{[6]} reproducing kernel particle method,^{[7–9]} smoothed particle hydrodynamics,^{[10,11]} moving particle semi-implicit method,^{[12]} and corrective smoothed particle method.^{[13]}

Smoothed particle hydrodynamics (SPH)^{[14]} is a transient, Lagrangian method, proposed by Lucy,^{[15]} Gingold and Monaghan^{[16]} in 1977. The SPH was originally invented for modeling the astrophysical phenomena in three-dimensional (3D) open space, and now has been widely applied to fluid dynamics such as numerical simulation of quasi-incompressible flows,^{[17]} multi-phase flows^{[18]} and heat conduction.^{[19]} In the SPH, a set of particles associated with physical properties are used to represent the state of the system. Since these particles can move within the computational domain in the Lagrangian frame, SPH is an attractive tool to model complex flows with free surfaces, moving interfaces and large deformations which are difficult for grid-based methods. Though the SPH has a very wide range of applications in incompressible flows,^{[20–22]} applications in heat transfer problems are still not mature enough.

Natural convection in a closed square cavity is a classic convection heat transfer problem, which is important not only for understanding heat transfer phenomena but also for studying problems such as cooling of electrical equipment, nuclear reactor systems and crystal growth in liquids.^{[23]} Hence, studies of heat transfer problems and obtaining their accurate solutions are still significant. In 1998, by using the SPH, Cleary^{[24]} modeled natural convective flows in a differentially heated square cavity for low and moderate Rayleigh numbers (10^{2}, 10^{4}, 10^{5}). The results showed that the SPH is sufficient to give excellent isotherms for low *Ra*, but when *Ra* = 10^{5}, obvious deviations appear between SPH results and FEM results. Chaniotis *et al*.^{[25]} used the remeshed SPH to simulate natural convection in a differentially heated cavity, displayed the contours of temperature and vorticity for *Ra* = 10^{3} and 10^{4}, and the relative error between the remeshed SPH solution and the benchmark solution of de Vahl Davis is less than 8%. Szewc *et al*.^{[26]} modeled natural convection with SPH based on non-Boussinesq formulation, the local Nusselt number distributions look very realistic, but contours for temperature and velocity are not smooth. Though it has been widely used in the field of solid mechanics and fluid mechanics, the conventional SPH is still poor in stability and accuracy.

One of the factors causing SPH instability is the non-uniform particle distribution.^{[27]} In the Lagrangian numerical simulation, particles move along the streamlines, which will cluster or scatter under certain conditions, such as natural convection. What is more, the non-physical behavior of the kernel function will weaken the interaction between particles when particle distances are within a close range. Then particles will continue to cluster, which will reduce the precision or interrupt the numerical simulation, so it is necessary to correct the particle position. In 2008, Nester *et al*.^{[28]} proposed a particle velocity correction, and applied this method to the finite volume particle method (FVPM) to maintain the uniformity of the moving particle cloud. In 2009, particle velocity correction was improved by Xu *et al*.,^{[29]} and the new approach called shifting particle position maintains accuracy and stability when simulating Taylor–Green vortices and vortex spin-down problems. In 2012, Hashemi *et al*.^{[30]} applied shifting particle position to SPH, the numerical results demonstrate that shifting particle position can alleviate defects produced by non-uniform distribution of particles and improve the stability of the numerical method. Shifting particle position not only needs less computational time but also effectively maintains the uniformity of particle distribution. Hence, this paper uses shifting particle position to weaken the non-uniform particle distribution effects on numerical simulation.

At present, in order to overcome the low accuracy of the conventional SPH, a large number of methods have been proposed. Finite particle method (FPM)^{[31]} is an enhanced version of SPH, in which the function and its first derivatives are calculated simultaneously to restore the particle consistency and improve the approximation accuracy. Zhang and Batra^{[32]} used numerical tests to show the FPM superiority over the CSPM, and the simulation results of the 2D transient heat conduction problems demonstrate that FPM is more accurate than CSPM. Liu *et al*.^{[31]} modeled the classic Poiseuille flow, Couette flow, shear driven cavity and a dam collapsing problem by FPM, and the numerical examples showed that the FPM is a very attractive alternative for simulating incompressible flows. On the basis of FPM, Huang *et al*.^{[33]} substituted the product of particle distance and kernel function for kernel gradient in FPM since the kernel gradient is not necessary in the whole computation, and this new approach is referred to as a kernel gradient free-SPH (KGF-SPH). KGF-SPH widens the selection of the kernel function, improves the simulation precision and simplifies the solving process because of the symmetric property of the corrective matrix. However, the conventional FPM and KGF-SPH differentiate the first derivatives to obtain the second derivatives, and the errors of the first derivatives will be brought into the second derivatives, which will reduce the accuracy of the simulation results.

In the present paper, the Improved KGF-SPH is proposed, which modifies the two-pass model in KGF-SPH with a new kind of discretization scheme for the Laplacian. Simple (1D and 2D heat conduction) and complex (natural convection in a closed square cavity) heat transfer problems are simulated by the Improved KGF-SPH. The rest of this paper is organized as follows. The conventional SPH, KGF-SPH, and Improved KGF-SPH are introduced in Section 2. In Section 3, the Improved KGF-SPH formulations for the governing equations and the state equation are briefly depicted. In Section 4, the particle shifting position is presented in detail. In Section 5, the 1D and 2D heat conduction problems are adopted to test the precision and stability of the Improved KGF-SPH. In Section 6, natural convection in a closed square cavity is modeled by the Improved KGF-SPH with shifting particle position. In Section 7, some conclusions are drawn from the present research results.

The concept behind SPH is based on an interpolation scheme, and in SPH, any function or physical property can be approximated by its value at a set of disordered points. The formulations of SPH are often divided into two key steps: the first step is the integral representation or the so-called kernel approximation of physical properties, the second one is the particle approximation of physical properties. In the first key step, for an arbitrary function *f*, the kernel approximation of the function value *f* (* r _{i}*) and its first derivatives ∇

*f*(

*) at*r

_{i}*=*r

*are respectively*r

_{i}^{[14]}

*f*(

*) is the function value at*r

_{j}*=*r

*,*r

_{j}*W*(

*−*r

_{i}*,*r

_{j}*h*) is the kernel function,

*h*is the smoothing length, ∇

_{i}

*W*(

*−*r

_{i}*,*r

_{j}*h*) is the kernel gradients, and d

*V*is the volume element.

In the particle approximation, the computational domain is discretized by a series of particles with some attributes such as mass, position, momentum and energy. Equations (

*f*=

_{i}*f*(

*),*r

_{i}*f*=

_{j}*f*(

*),*r

_{j}*W*=

*W*(

*−*r

_{i}*,*r

_{j}*h*), (∇

*f*)

_{i}= ∇

*f*(

*), and ∇*r

_{i}_{i}

*W*= ∇

_{i}

*W*(

*−*r

_{i}*,*r

_{j}*h*). The subscript

*j*represents the neighboring particle of particle

*i*,

*m*and

_{j}*ρ*

_{j}represent mass and density of particle

*j*respectively.

In 2D case, the most widely used SPH discretization scheme for the Laplacian is^{[34]}

*α*and

*β*are the dimension indices repeated from

*x*to

*y*,

*r*and

^{α}*r*represent the components of

^{β}*r*

*α*and

*β*directions respectively,

*δ*is the Dirac delta function: when

^{αβ}*α*=

*β*,

*δ*= 1, when

*α*≠

*β*,

*δ*= 0. According to Eq. (

In 2015, based on FPM, Huang *et al*.^{[33]} developed a kernel gradient free-SPH (KGF-SPH) which only involves kernel function itself in kernel and particle approximation. In the KGF-SPH, performing Taylor series expansion at a certain point * r* =

*and retaining the first order derivatives, an arbitrary function*r

_{i}*f*(

*) can be expressed as follows:*r

_{j}*o*(

*h*

^{2}) is the remainder of the expansion. Multiplying both sides of Eq. (

*W*and (

*−*r

_{j}*)*r

_{i}*W*respectively, ignoring the remainder of the expansion and integrating over the computational domain can yield the following equations:

In 2D space, equation (

*A*

*x*=

_{ji}*x*

_{j}−

*x*

_{i}and

*y*=

_{ji}*y*

_{j}−

*y*

_{i}.

It can be observed from Eq. (

The two-pass model, which obtains second derivatives by nesting two first derivative operations, is used to discretize the Laplacian in the conventional KGF-SPH. So in 2D space, the KGF-SPH discretization scheme for the Laplacian is

*f*

_{j,x}and

*f*

_{j,y}are the first derivatives of particle

*j*in the

*x*and

*y*directions respectively. The subscript

*k*represents the neighboring particle of particle

*j*,

*x*

_{k j}=

*x*

_{k}−

*x*

_{j},

*y*=

_{k j}*y*−

_{k}*y*.

_{j}In KGF-SPH, the Laplacian is calculated by using Eq. (

In order to overcome these two drawbacks, in this paper we propose a new discretization scheme to approximate the Laplacian, and apply it to the KGF-SPH. The following section provides the derivation process of the new Laplacian discretization scheme in detail.

In 2D space, performing Taylor series expansion at a nearby point (*x _{i}*,

*y*

_{i}) and retaining the second derivatives, a function

*f*at point (

*x*,

_{j}*y*

_{j}) can be expressed as follows:

*f*

_{i,xx}and

*f*

_{i,yy}are the second derivatives of particle

*i*in the

*x*and

*y*directions respectively,

*f*

_{i,xy}is the second order mixed partial derivative of particle

*i*.

Multiplying both sides of Eq. (*W*, and transposing can yield the following equation:

Equations (

For natural convection in a closed square cavity, assuming the dynamic viscosity coefficient *μ* to be a constant, the governing equations in the Lagrangian frame are^{[14]}

*t*represents the time,

*ρ*is the density,

*x*and

*y*are the spatial coordinates,

*u*and

*v*are the velocity components in the

*x*and

*y*directions respectively,

*p*is the pressure,

*T*represents the temperature, thermal diffusivity

*α*=

*λ*/(

*ρc*),

_{p}*λ*being the thermal conductivity, and

*c*

_{p}the specific heat capacity.

Using Eqs. (

*u*=

_{ji}*u*−

_{j}*u*

_{i},

*v*=

_{ji}*v*−

_{j}*v*

_{i}, and

*p*=

_{ji}*p*

_{j}−

*p*.

_{i}The above governing equation is not closed, so the state equation should be added. In this paper, the following weakly compressible state equation is used^{[17]}

*ρ*

_{∞}is the reference density,

*c*is the artificial sound speed, and its value is ten times as large as the maximum velocity.

In order to enhance numerical stability, shifting particle position can be used to maintain the uniformity of particles. Shifting particle position employed in this paper modifies the particle position, and the position shift *δ r _{i}* reads

*C*is a constant: the bigger the

*C*value, the greater the effect of the method is, and

*C*is set to be 0.5 in this paper;

*α*is the shifting magnitude;

*R*

_{i}is the shifting vector of particle

*i*;

*α*and

*R*

_{i}are solved by the following equations:

*U*

_{max}is the maximum particle velocity, d

*t*is the time step,

*M*is the number of neighboring particles around particle

_{i}*i*;

*r*is the distance between particle

_{i j}*i*and particle

*j*;

*n*

_{ij}is the unit distance vector between particle

*i*and particle

*j*;

*r̄*is the average particle spacing in the neighborhood of

_{i}*i*, which reads

Since the particles are shifted, the physical properties carried by particles are also changed. The physical properties are corrected by the following equation:

*φ*is a general variable representing any physical properties,

*i*and

*i*′ are the particle old position and its new position respectively. Shifting particle position can obtain an isotropic particle distribution, so it is used to improve numerical stability of natural convection in a closed square cavity.

In order to validate the Improved KGF-SPH, two simple conduction heat transfer problems are studied in this section. The Improved KGF-SPH results are compared with those from SPH, KGF-SPH, and analytical solutions to show the high accuracy and stability of the Improved KGF-SPH.

The computational domain for 1D heat conduction in this paper is the line segment from 0 to 1, at *t* moment, the temperature *T*(*x*, *t*) is governed by the following heat conduction equation:

*T*is the temperature,

*t*is the time, and

*x*is the particle position. The initial and boundary conditions are

*T*(

*x*, 0) = sin

*πx*and

*T*(0,

*t*) =

*T*(1,

*t*) = 0, respectively. In this case, the analytical solution is given by

Figure *Δ* = *L*/50 = 0.02, the smoothing length is *h* = 1.3*Δ*. According to the CFL condition, the time step is set to be d*t* = 2 × 10^{−4} s.

Figure *t* = 0.2 s for 1D heat conduction, computed by SPH, KGF-SPH, and Improved KGF-SPH with the analytical solution. It can be seen from Fig.

In this section, the Improved KGF-SPH formulations are applied to the 2D heat conduction case with a simple square cavity geometry whose side length is *L*. At *t* moment, the temperature *T*(*x*,*y*,*t*) is governed by

*T*(

*x*,

*y*,0) = 0, and

*T*(0,

*y*,

*t*) =

*T*(0.1,

*y*,

*t*) =

*T*(

*x*,0,

*t*) =

*T*(

*x*,0.1,

*t*) = 1, respectively.

Figure *L* = 0.1 m, fluid particles are uniformly distributed in the computational domain, and virtual particles are distributed evenly on the boundary. The smoothing length is *h* = 1.3*Δ*, *Δ* is the particle spacing, and the time-step is d*t* = 10^{−6} s. As shown in Fig. *A*(0.005,0.01) is the monitor particle, so in this paper we will monitor the temperature history for point *A*.

In order to verify our numerical method, the solution of the 2D heat conduction problem by the finite particle method (FVM) with 101×101 nodes is taken as a reference solution. For the SPH, KGF-SPH, and Improved KGF-SPH, 11×11 and 21×21 uniformly distributed particles are employed.

Figure *x* = 0.05 m at three times, i.e., *t* = 1.5 × 10^{−4}, 3.0 × 10^{−4}, and 4.5 × 10^{−4} s, computed with SPH, KGF-SPH, Improved KGF-SPH, and FVM. Because of the symmetry of the problem about the lines *x* = 0.05 m and *y* = 0.05 m, temperature distribution in the region *y* ∈ [0.05,0.1 m] is depicted. Figure

In order to show the Improved KGF-SPH superiority over SPH and KGF-SPH clearl*y*, figure *y*, the results of KGF-SPH will not be exhibited in the following sections.

The Improved KGF-SPH with shifting particle positions is used to model natural convections in a closed square cavity at different Rayleigh numbers in this section. In order to verify our Improved KGF-SPH formulation, the Improved KGF-SPH results are compared with the corresponding SPH and FVM simulation results.

Schematic diagram for natural convection in a closed square cavity is provided in Fig. *L*, both the upper and lower horizontal walls are adiabatic, the left and right vertical walls are isothermal, maintained respectively at the maximum temperature *T*_{H} and the minimum temperature *T*_{C}. Throughout the simulations, gravity is in negative *y* direction and Boussinesq approximation is used. The Boussinesq approximation means that the density differences are confined to the buoyancy term *gβ*_{T}(*T* − *T*_{0}), where *g* represents the gravitational acceleration, *β*_{T} is the coefficient of thermal expansion, and the reference temperature *T*_{0} = (*T*_{H} + *T*_{C})/2.

Particle distribution for natural convection in a closed square cavity is shown in Fig. *h* = 1.3*Δ*, where *Δ* is the particle spacing, and the time step is d*t* = 5 × 10^{−6} s.

Natural convection in a closed square cavity is studied at *Ra* = 10^{4}, 10^{5} and 10^{6} to validate the Improved KGF-SPH. In this case, the Prandtl number *Pr* = *μ*/(*ρ α*) = 0.71, the dynamic viscosity coefficient *μ* = 1.707 × 10^{−5} kg/m·s, density *ρ* = 1.225 kg/m^{3}, the thermal diffusion coefficient *α* = 1.963× 10^{−5} m^{2}/s, the thermal expansion coefficient *β*_{T} = 4.095 × 10^{−3} K^{−1}. *T*_{H} = 283 K, *T*_{C} = 273 K, and the initial temperature of the fluid inside the square cavity is *T*_{0}. The reference density *ρ*_{∞} in state equation is set to be 1.225 kg/m^{3}, and for convenience, the artificial sound speed *c* values are set to be 0.25 m/s, 0.4 m/s and 0.8 m/s for *Ra* = 10^{4}, 10^{5}, and 10^{6} respectively. In this paper, the size of Rayleigh number Ra is adjusted by changing the sizes of *g* and *L*. In the numerical simulation of natural convection in a closed square cavit*y*, shifting particle position is adopted to improve the numerical stability.

There are two kinds of boundary conditions in natural convection problems: adiabatic boundaries and isothermal boundaries. In this paper, we use the projecting method proposed by Huang *et al*.^{[33]} to calculate physical properties of virtual particles. The projecting method first finds out the normal direction *n*_{vp} of all virtual particles’ position vector * r _{vp}*, and extends virtual particles by an increment

*η Δ*into the flow field along the normal direction

*n*

_{vp}. Then we can obtain new points located on

*, which are the projecting points corresponding to the virtual particles.*r

_{pp}*η*is an adjustable parameter,

*Δ*is the initial particle spacing, the subscripts

*vp*and

*pp*represent virtual particles and projecting points respectively. In order to understand the projecting method, the value of

*η*is given through four specified virtual particles

*A*,

*B*,

*C*, and

*D*on each side as shown in Fig.

The adiabatic boundary condition is a Neumann-type boundary condition with *∂T*/*∂ n* = 0, where

*is the normal direction of the wall. On the adiabatic boundaries, the position vectors of virtual particles and projecting points satisfy the following relations:*n

*x*and

*y*represent the components of the position vector. The principle of selecting

*η*value is that the distances of projecting points to the walls are all

*Δ*/2. For example, as shown in Fig.

*A*is located on point

*M*, and the corresponding projecting point of particle

*B*is located on point

*N*, thus the value of

*η*is −1.5 and 3.5 respectively. After finding the projecting points, the physical properties of virtual particles are calculated by the following equations:

*j*is the neighboring particle of projecting point

*pp*, only fluid particles are involved. The pressure of virtual particles

*p*

_{vp}is calculated by Eq. (

The isothermal boundary condition is a Dirichlet-type boundary condition. On the isothermal boundaries, the position vectors of virtual particles and projecting points satisfy the following relations:

*η*value is that the distances of projecting points to the walls are all

*Δ*/2. For example, as shown in Fig.

*C*is located on point

*E*, and the corresponding projecting point of particle

*D*is located on point

*F*, thus the value of

*η*is 4.5 and –2.5 respectively. The physical properties of virtual particles are solved by the following equations:

*j*is the neighboring particle of projecting point

*pp*, only fluid particles are involved,

*T*

_{W}is the temperature of the walls,

*T*

_{W}=

*T*

_{H}on the left vertical wall,

*T*

_{W}=

*T*

_{C}on the right vertical wall. The pressure of virtual particles

*p*

_{vp}is calculated by Eq. (

For natural convection in a closed square cavity with *Ra* = 10^{4}, *g* = 8.350 m/s^{2}, *L* = 0.02 m, 149×149 fluid particles are evenly placed in the computational domain. Figures *t* = 4.5 s, obtained by (a) FVM, (b) SPH, and (c) Improved KGF-SPH, respectively. In Fig. *Ra* = 10^{4}.

The simulations for *Ra* = 10^{5} are performed with 149×149 fluid particles, *g* = 5.344 m/s^{2}, *L* = 0.05 m, and the temperature contours t *t* = 20 s are presented in Fig. *u* contours and vertical velocity *v* contours obtained using (a) FVM, (b) SPH and (c) Improved KGF-SPH at *t* = 20 s for the natural convection in a closed square cavity with *Ra* = 10^{5}. Though temperature contours obtained by SPH accord well with that of FVM, the velocity contours obtained by SPH are poor in stability. As shown in Figs.

For *Ra* = 10^{6}, *g* = 6.680 m/s^{2}, *L* = 0.1 m, the boundary layer becomes thinner with Rayleigh number increasing, so more particles (199×199 fluid particles) are used to model the natural convection in a closed square cavity. Figure *t* = 30 s, obtained by using (a) FVM, (b) SPH, and (c) Improved KGF-SPH. Both the SPH and Improved KGF-SPH results accord well with that of FVM. The corresponding horizontal velocity *u* and vertical velocity *v* contours are depicted in Figs. *y*, but the Improved KGF-SPH results are still very smooth and in good agreement with those obtained by FVM.

There are two distinct flow patterns for natural convection in a closed square cavity: heat conduction and heat convection. It can be seen from Figs. *Ra* increases, the interior temperature contours become more and more horizontal, since the heat convection dominates the higher-*Ra* range. Nevertheless, the contours near the adiabatic boundaries are normal to the horizontal walls. In Figs. *Ra*, these two eddies move closer to the adiabatic walls. Similarl*y*, two vertical eddies are visible in the vertical velocity contours (Figs.

It can be seen from Fig. *Ra* number increases, the boundary layer becomes thinner, and the flow inside the cavity becomes more and more complex. The temperature contours obtained by using the SPH and Improved KGF-SPH are in good agreement with that of FVM, so both the SPH and the Improved KGF-SPH can simulate the temperature field for the natural convection exactly. However, the velocity contours obtained by using the SPH become poorer for higher Ra, especially when *Ra* = 10^{6}, but the Improved KGF-SPH results are in good agreement with FVM results. To sum up, the stability of the Improved KGF-SPH is better than SPH.

In order to validate the accuracy of Improved KGF-SPH, the velocity profiles obtained by SPH and Improved KGF-SPH are compared with those of FVM. For convenience, the dimensional variables are non-dimensionalized by the following equations: *x** = *x*/*L*, *y** = *y*/*L*, *u** = *uL*/*α*, *v** = *vL*/*α*. The non-dimensional horizontal velocity *u** distribution along the mid-width (*x* = *L*/2) and non-dimensional vertical velocity *v** along mid-height (*y* = *L*/2) over the *Ra* range of 10^{4} to 10^{6} are presented in Figs. *Ra* = 10^{6} is also presented in Fig. *Ra* = 10^{4} and 10^{5}, but when *Ra* = 10^{6}, there is a large error between SPH and FVM results while the Improved KGF-SPH results are still close to the FVM results. In Fig. *Ra* = 10^{6}, we can see that the maximum non-dimensional vertical velocity obtained using Improved KGF-SPH is closer to FVM results than that of SPH. In general, Improved KGF-SPH is more accurate than SPH for modeling natural convection in a closed square cavity.

The maximum non-dimensional horizontal velocity ^{[35]} and Wan *et al*.^{[36]} are also summarized in Table *Ra* = 10^{6} with error about 2.66%, so the Improved KGF-SPH is more accurate than SPH in modeling natural convection in a closed square cavity.

A new Laplacian discretization scheme is proposed to avoid bringing errors of the first derivatives to the second derivatives in the KGF-SPH kernel and particle approximation. The superiority of the new Laplacian discretization scheme applied in KGF-SPH is verified by modeling 1D and 2D heat conduction problems. The Improved KGF-SPH is successfully used to study natural convection in a closed square cavity with *Ra* = 10^{4}, 10^{5}, and 10^{6}. It can be seen that the SPH overestimates the heat transfer, the KGF-SPH results oscillate violentl*y*, while the Improved KGF-SPH results accord well with the analytical solution and the reference solution. The Improved KGF-SPH is more accurate than SPH and stabler than KGF-SPH, thus it is an attractive tool to study and model the heat transfer problems. In this paper, only 1D and 2D simulations are done; as a next step, our method can be extended to solving 3D heat transfer problems. Furthermore, the application of the proposed method can be broadened to other problems controlled by governing equations, such as viscous fluid flow.

**Reference**