Heat transfer is a fundamental physical process that plays a critical role in a wide range of engineering applications, from manufacturing and materials processing to energy systems and modern thermal technologies. These include heating plate panels, melting and freezing, metal casting, and thermal spraying [1, 2, 3, 4, 5]. Understanding the mechanisms of heat transfer is essential to improve energy eficiency, predict the thermal response of a system, and design devices that rely on temperature control.

Various experimental studies have been conducted to observe the conduction patterns and temperature distribution in solid materials [6, 7, 8]. However, experimental approaches are generally expensive and time-consuming. Therefore, mathematical models have become important alternatives for systematically studying heat transfer phenomena. In this study, the heat transfer process is modeled using the heat equation, which represents the most fundamental form of difusion [9]. Conceptually, the heat equation belongs to the family of difusion-type partial diferential equations (PDEs), including the Extended Fisher–Kolmogorov equation (EFK) [10, 11] and the Kuramoto–Sivashinsky equation (KS) [12, 13], which are fourth-order equations.

The analytical solution of the heat equation depends heavily on the boundary conditions and becomes intractable when the geometry or boundaries are sufficiently complex. In such cases, numerical approaches provide a more practical framework for modeling realistic heat transfer scenarios [14, 15]. Classical gridbased schemes, such as the Finite Diference Method (FDM) [16], the Finite Volume Method (FVM) [17], and the Finite Element Method (FEM) [18], remain widely used for their robustness in structured domains. Various advanced discretization techniques have also been developed, including the One-Step One-Hybrid Block Method for nonlinear oscillatory PDEs [19], Hermite spline methods for the KS equation [20], and hybrid quintic Hermite spline collocation for the EFK equation [21]. Despite their strengths, grid-based methods can become restrictive when dealing with highly irregular boundaries, large deformations, or free-surface problems, motivating the use of more flexible mesh-free particle approaches such as Smoothed Particle Hydrodynamics (SPH).

SPH was first introduced by Lucy [22] and later developed independently by Gingold and Monaghan [23] for astrophysical simulations. Since then, it has been widely used to simulate solid mechanics and fluid dynamics [24, 25, 26]. Unlike grid-based schemes that require a predefined mesh, the SPH represents the domain as a set of interacting particles that carry quantities such as mass, position, velocity, and temperature. This mesh-free nature allows SPH to naturally handle large deformations, moving boundaries, and highly irregular geometries. The governing PDEs are reformulated as ordinary diferential equations using integral approximations based on kernels, which form the core of the SPH methodology [27].

The use of SPH for heat transfer problems began with the pioneering work of Cleary and Monaghan [28], who demonstrated its ability to solve heat conduction equations. Since then, SPH has been extensively applied to a broad range of thermal processes, including conduction, convection, and radiation in complex geometries [29, 30, 31, 32, 33, 34, 35]. Recent developments have further expanded its applicability to enhanced heat transfer in nanofluids [36] and chemically reactive systems where thermal difusion is coupled with reaction kinetics [37, 38]. Recent advances in manufacturing further emphasize the importance of modeling internal heat sources, as demonstrated by Lin et al. [39]), who applied a ray-tracing SPH model to represent laser-induced heating. In most classical heat transfer studies, the heat source is imposed directly within the heat equation [40, 41].

Although SPH has been applied in various heat transfer studies, several research gaps remain. The comprehensive validation of SPH against two-dimensional analytical solutions, especially for Dirichlet and Neumann boundary conditions, is still limited. In addition, the efect of the geometry of the internal heat source on the heat transfer rate in isolated plates, which is analytically challenging, has not been widely examined using SPH. To address these gaps, this study provides a systematic validation of the SPH against two-dimensional analytical solutions under Dirichlet and Neumann boundary conditions. Furthermore, we applied SPH to simulate heat transfer in plates containing internal heat sources of diferent geometries but equal area, allowing the assessment of how geometric variations influence heat transfer rates. The standard heat equation without a source term is employed, with the internal heat source imposed through Dirichlet conditions within the source region.

The remainder of this paper is organized as follows. Section 2 presents the mathematical formulation of the two-dimensional heat equation. Section 3 describes the SPH numerical scheme used to approximate the solution. Section 4 provides model validation by comparing the SPH results with analytical solutions. Section 5 investigates the influence of the geometry of the internal heat source on the heat transfer rate. Finally, Section 6 summarizes the main findings of this study.

This section discusses the governing equations for heat transfer in a plate. Consider a thin plate made of heat-conducting material. Let be the temperature of the plate at point and time . This study is restricted to a rectangular plate , where .

Under the ideal assumption of uniform density and no external heat sources, u satisfies the two-dimensional heat equation:

for . Thermal difusivity α is defined as

where is the thermal conductivity, is the heat capacity, and is the density of the material. In this study, the physical domain consists of a 10 cm × 10 cm carbon steel plate characterized by the material properties , and . In steady state, equation (1) reduces to Laplace’s equation.

The analytical solution of the two-dimensional heat equation (1) with generalized Dirichlet boundary conditions has been discussed by Hsu et al. [42]. As for Neumann boundary conditions, to the best of our knowledge, the analytical solutions have not been discussed in the literature. However, Subani et al. [43] have addressed this for the one-dimensional case, where the heat equation (1) in one-dimensional form with homogeneous Neumann boundary conditions is solved in two stages: variable separation and the superposition principle.

In the following, we present the analytical solution of the two-dimensional heat equation (1) with homogeneous Neumann boundary conditions. The plate temperature at the beginning of the observation is given as the initial condition:

Furthermore, the homogeneous Neumann boundary conditions are as follows:

Physically, this means isolating the four sides of the plate so that no heat enters or leaves the plate through these sides.

Using a procedure similar to [43], we obtain an analytical solution of the heat equation (1) with the initial conditions (3) and the boundary conditions (4) as follows:

where

The detailed derivation of this analytical solution is provided in the Appendix.

In this section, we describe the SPH formulation used to approximate the twodimensional heat equation (1). The domain is discretized into a number of particles with the initial positions of the -particles being . Two types of initial particle arrangement are commonly used: structured and unstructured, as illustrated in Figure 1. A structured layout produces higher operator accuracy [44], but is less flexible to capture intricate geometries. In contrast, an unstructured distribution allows particles to reach narrow features of the domain, although the resulting particle disorder may introduce extra numerical noise [45].

The function at and time can be approximated by volume integration using the kernel function (also called the weight function) , as follows

for , where is the distance between the particles and , and are the smoothing lengths. Consistent with typical SPH settings [24, 27, 45], we define the smoothing length as , where is the initial spacing of the particles.

Neighborhood particles are a crucial component in particle kernel approximation. Define the vector as the distance vector between particles and ,

The distance between the target particle and the neighboring particle is defined using the Euclidean norm

In this study, we use a multiplier on the smoothing length. If , then the particle is considered a neighbor of the particle . For example, Figure 2 (left) illustrates the definition of neighboring particles, where particle has 20 neighboring particles.

Furthermore, several smoothing kernels are commonly used, including cubic spline (B-spline), Gaussian, Supergaussian, Quadratic, and Quintic kernels. The kernel function used in this study is the cubic spline, which is given as

where . The constant varies for each dimension to satisfy the kernel normalization requirements. For a two-dimensional system, . Note that the kernel (7) is a smooth polynomial function, which allows it to be diferentiated at all . The cubic spline kernel function and its derivatives are plotted in Figure 2 (right). It is also important to note that the values of the kernel function and its derivatives are zero for 2 \end{document} ]]>, meaning that particles located more than away from the target particle do not contribute to the approximation.

Since the domain is modeled as discrete particles, the continuous integral (6) needs to be converted to a discrete sum over all particles in the domain. Let be the set of neighboring particles surrounding the target particle , which is discretized as

The approximation (6) can be expressed in its discrete particle form as follows:

where and are the mass and density of the particle , respectively.

The kernel function derivative is used for the Laplace approximation. The kernel gradient vector is given as

where

Finally, the gradient and Laplace approximations of the SPH method are given as follows:

After approximating the Laplace form using the SPH method, we proceed with the time integration of the heat equation (1). To ensure a stable solution, the time step must satisfy the Courant-Friedrichs-Lewy (CFL) condition. This condition ensures that the numerical domain of dependence contains the physical domain of dependence, which means that the numerical propagation speed must not be smaller than the physical propagation speed [46, 47]. In the case of twodimensional heat difusion, the time step is formulated as

The temperature of the particle at time is then computed using the forward Euler method, given by

To handle the boundary value problem, SPH utilizes the concept of ghost particles, which are additional particles located outside the domain or around the boundary; refer back to Figure 1. These particles are not included in the SPH computation, but are arranged to approximate the prescribed boundary conditions. Dirichlet boundary conditions are enforced by setting the temperature of the ghost particles equal to the specified boundary values, as described in [32]. Implementing homogeneous Neumann boundary conditions poses additional challenges, such as assigning the temperatures of ghost particles so that the normal temperature gradient vanishes, which requires arranging the particles near the boundary in a structured manner so that the ghost particle temperature can be mirrored from the interior particle across the boundary.

In summary, the steps for using the SPH method to simulate heat transfer are outlined in Algorithm 1.

All simulations were performed on a computer equipped with an AMD Ryzen 9 PRO 5945 processor (12 cores, 24 CPUs, 3.0 GHz), 128 GB RAM, and a Windows 11 operating system. Numerical calculations were performed using MATLAB R2024b.

In this section, the SPH results obtained using the numerical procedure described in Algorithm 1 are compared with the corresponding analytical solutions to assess the accuracy of the method. The computational domain is expressed as . The plate is discretized using a structured particle arrangement, where and denote the initial particle spacing in the and directions, respectively, and are taken to be equal in this study.

The agreement between the SPH and the analytical reference is quantified using the normalized root mean square error (nRMSE). This metric is chosen because it provides a scale-independent measure of accuracy, allowing the error to be interpreted relative to the thermal range of the system and enabling fair comparison between diferent test cases. The nRMSE is defined as

where and denote the maximum and minimum values of the analytical solution in the domain, respectively.

Finally, SPH is used to simulate heat transfer in more complex cases. We investigate the efect of an internal heat source located in the center of the plate on the rate of heat spreading. Heat sources with various geometries, all having the same area, are considered. The objective is to analyze which geometric shapes distribute heat more rapidly than others. All four sides of the plate are insulated, which prevents heat from entering or exiting the edges. This case is quite complex to solve analytically, as it involves nonhomogeneous Dirichlet and homogeneous Neumann boundary conditions, and the geometric shape of the heat source is not simple. To the best of our knowledge, this issue has not been previously discussed in the literature.

We tested six heat source geometries: circular, triangular, square, hexagonal, star, and star of David. The radius of the circular heat source is 1 cm, and all heat sources have the same area. The dimensions of these heat sources are shown in Figure 8. Initially, the plate has a uniform temperature of . Then a constant heat source of is applied.

To allow particles to reach narrow corners, we use an unstructured particle distribution, and the particle spacing is set to be no larger than 0.2 cm. The simulation is carried out up to using a time step of . The number of particles and the corresponding computational time for each heat-source eometr are summarized in Table 2.

The heat distribution at for all heat sources is shown in Figure 9. It can be seen that the heat spreads radially toward the edges of the plate, forming a circular pattern. The temperature decreases monotonically from the center to the edges. It can be observed that, among all heat sources, those with star geometries disperse heat more rapidly.

To compare the rate of heat propagation, we plot the contours of a selected temperature for all heat source geometries. Here, we plot the contour for , which is depicted in Figure 10 (left). It can be seen that the temperature range for star geometries is wider than for other geometries. In addition, we also calculated the average plate temperature outside the heat source , which is shown in Figure 10 (right). The average temperature also shows that the star geometry disperses heat more quickly, compared to the other five geometries. The star geometry spreads heat faster because its extended tips reach farther toward cooler areas of the plate, creating stronger temperature diferences that accelerate difusion.

This study demonstrates that the Smoothed Particle Hydrodynamics (SPH) method can accurately simulate two-dimensional heat transfer on thin plates under various boundary conditions. Validation against analytical solutions for Dirichlet and homogeneous Neumann problems shows excellent agreement, as indicated by consistently low nRMSE values. The method was also applied to plates containing internal heat sources of equal area but diferent geometries, revealing that starshaped sources promote faster heat spread than other configurations. These results highlight the flexibility of SPH in handling complex boundary specifications and geometric variations.

Despite these strengths, the present formulation has several limitations. SPH requires user-defined parameters such as the smoothing length and the kernel function, both of which introduce uncertainty. In addition, the model used here solves only a scalar temperature field with fixed particle positions, which does not exploit the full potential of SPH for deforming or moving domains. The governing physics is limited to pure conduction, and the validation is restricted to rectangular domains and analytical benchmarks. Future work may include validation against analytical solutions in circular domains as well as comparisons with laboratory experiments, which will require more complex governing equations that incorporate heat sources, convection, or phase change. Extensions to fully coupled thermofluid systems, moving particles, and three-dimensional geometries would further strengthen SPH as a versatile tool for realistic heat-transfer modeling.

Details of the search for analytical solutions of the two-dimensional heat equation for homogeneous Neumann boundary conditions: