Second-Order Symmetric Smoothed Particle Hydrodynamics Method for Transient Heat Conduction Problems with Initial Discontinuity

Abstract: To eliminate the numerical oscillations appearing in the first-order symmetric smoothed particle hydrodynamics (FO-SSPH) method for simulating transient heat conduction problems with discontinuous initial distribution, this paper presents a second-order symmetric smoothed particle hydrodynamics (SO-SSPH) method. Numerical properties of both SO-SSPH and FO-SSPH are analyzed, including truncation error, numerical accuracy, convergence rate, and stability. Experimental results show that for transient heat conduction with initial smooth distribution, both FO-SSPH and SO-SSPH can achieve second-order convergence, which is consistent with the theoretical analysis. However, for oneand two-dimensional conduction with initial discontinuity, the FO-SSPH method suffers from serious unphysical oscillations, which do not disappear over time, and hence it only achieves first-order convergence; while the present SO-SSPH method can avoid unphysical oscillations and has second-order convergence rate. Therefore, the SO-SSPH method is a feasible tool for solving transient heat conduction problems with both smooth and discontinuous distributions, and it is easy to be extended to high dimensional cases.


Introduction
Heat transfer is generally defined as the phenomenon where heat is spread from one system or one part to another, and can be divided into three modes: heat conduction, heat convection, and thermal radiation.As the main form of heat transfer in solids, heat conduction is applied in many industrial fields, such as rubber vulcanization, steel forgings heating, and the design of thermal conductors and insulation equipment.It is of high importance to industrial production to explore the inherent law of heat conduction and to analyze the temperature distribution in heat conduction processes.
Theoretical methods cannot, in general, solve these complex heat transfer problems, while physical experiments suffer from high measurement cost.Based on computational simulations the Numerical Heat Transfer (NHT) technology [1] has been rapidly developed and the related theories have made great progress.Numerical methods for solving heat transfer problems are mainly mesh-based methods, such as finite difference methods (FDM) [2], finite volume methods (FVM), and finite element methods (FEM) [3].In recent years, due to the inherent advantages, meshless methods have received much attention and experienced rapid development.The major difference between mesh-based methods and meshless ones is that, in meshless methods, only node information is needed and the predefined node connectivity is not necessary.Zhang et al. [4] applied the Lattice Boltzmann method (LBM) to solve the nonlinear heat conduction equation.Singh et al. [5] used the Element-Free Galerkin (EFG) method to simulate transient heat conduction problems in isotropic homogeneous objects.As one of the earliest and most powerful meshless methods, the smoothed particle hydrodynamics (SPH) [6,7] method invented for astrophysical problems has been extended to complex fluid dynamics problems, especially for those with free surface and moving interfaces.Excellent reviews of the development of SPH and its diverse applications have been provided [6,[8][9][10][11][12][13][14].The development of SPH for heat conduction problems was mainly completed after the 2000s; we review the relative works below.It is worth noting that the SPH formulations for heat equations are by no means complete, and we will focus on the most representative and promising ones.Although great success and diverse applications have been achieved, there exist two inherent drawbacks [6] in the traditional SPH method: first, the boundary deficiency due to truncated kernel support of particles close to the boundary and second, accuracy degeneration when the particle distribution becomes highly irregular.The accuracy degeneration problem is closely related to the lack of completeness of the SPH approximations [15].To remedy these two drawbacks, several corrective techniques have been proposed, generating many corrective methods such as the corrective smoothed particle method (CSPM) [16], modified smoothed particle hydrodynamics (MSPH) [17] method, and moving least squares particle hydrodynamics (MLSPH) [18] method.
As conduction is often involved in other processes (e.g., viscous flow), there have been only a few attempts to entirely and exclusively solve the conduction problems using traditional SPH.For the discretization of the heat equation, there are three different ways.A simple discretization of the heat equation in SPH is obtained by firstly estimating the heat flux using the standard first derivative approximation, and then estimating the divergence in a second step.This method was originally used by Lucy [19] to simulate the thermal energy evaluation in the modeling of collapse and fission of optical thick protostar.The numerical results indicated that the log of density agreed to within approximately 1% of the exact solution.However, the disagreement between the exact solution and the SPH estimate of the divergence of the flux increased by 25% for particles close to the surface boundary.To remedy the boundary deficiency problem, Jeong et al. [20] deduced consistent estimation of near-boundary corrections for system variables using the ghost particle method and successfully verified this strategy against heat conduction tests.However, as shown by Brookshaw [21], this method was sensitive to distorted particle distribution, which can be traced to double differentiation of the kernel.Furthermore, another practical disadvantage is that, with an intermediate variable, the heat flux needs to be computed and stored in a separate SPH loop.The second method is the standard SPH kernel approximation of the second derivatives pioneered by Takeda et al. [22], which is sensitive to particle disorder [20] and may lead to divergence of the simulation.To remedy the divergence caused by distorted particle distributions, Chaniotis et al. [23] developed a periodic re-initialization (remeshing) technique.They found that the SPH methodology with remeshing procedure is capable of direct numerical simulation (DNS) quality simulations.In the above two methods, either a double calculation of the first derivative or a calculation of the second derivative of the kernel is needed.Both are sensitive to particle disorder that will most likely occur in fluid flows.It is hence advantageous to use a simpler SPH approximation of the Laplacian operator, which should only involve first-order derivatives of the kernel.This idea is commonly used now for heat conduction in SPH [10][11][12]14,[24][25][26][27][28][29].Such a method was pioneered by Brookshaw [21], further developed by Cleary [24,25], and its three-dimensional derivation was given by Jubelgas et al. [26] and Rook et al. [27].To take into account a discontinuous thermal conductively with continuous heat flux across the material interfaces, Cleary and Monaghan [25] proposed to replace the averaged conductivity by its harmonic mean, such that multiple materials with substantially different conductivities and specific heat can be accurately simulated [28,29].Based on Taylor series expansion both for regular and irregular particle distributions, Fatehi and Manzari [30] analyzed the truncation errors in the above three SPH approximations for second derivative and found that none of the three schemes has the first-order completeness.A new scheme for second derivatives was proposed using a modified renormalization tensor and applied to heat conduction problem.The key idea is similar as that in the MSPH method [17], where kernel derivatives are involved in the modifications.Recently, Francomano and Paliaga [31] highlighted the numerical insights in MSPH approximations and found that the request for greater accuracy needs kernel function derivatives with order up to the desired accuracy order in approximating the function or higher for the derivatives.
For heat conduction problems, there are few related studies using the corrective SPH methods.Early explorations are due to Chen et al. [16] and Zhang and Betra [17].Recently, Jiang et al. [32] presented the first-order symmetric smoothed particle hydrodynamics (FO-SSPH) method (named because the Taylor series of the function is expanded up to first derivatives), which has second-order accuracy and the symmetry property of the moment matrix.Furthermore, it does not involve kernel derivatives and can be directly applied to mixed boundary conditions.Although the FO-SSPH method can more effectively improve the numerical accuracy and stability than the standard SPH [32], it can only approximate the first-order spatial derivatives directly.Therefore, heat flux terms have to be introduced to transform the second-order heat equation into two first-order partial differential equations.Similar to the Forward Time Central Space (FTCS) method [33], FO-SSPH is only applicable for heat conduction problems with smooth temperature distribution.If the initial distribution has a discontinuity, unphysical oscillations will be inevitable and they cannot damp out naturally.This is because the temperature gradient at a discontinuous position tends to be infinite in theory, while the approximated first derivative in Taylor expansion can only be a finite value, so the introducing of the heat fluxes will magnify the error at this time.Therefore, it is difficult for FO-SSPH to solve transient heat conduction problems with initial discontinuity.Jiang et al. [32] concluded that the FO-SSPH has second-order convergence rate through numerical experiments.However, through theoretical analysis and numerical tests, we found that it is true for transient heat conduction problems with initial smooth distribution, but the FO-SSPH has only first-order convergence rate for problems with initial discontinuity.
Based on the above discussion, a second-order symmetric smoothed particle hydrodynamics (SO-SSPH) method (named because the Taylor series of the function is expanded up to second derivatives) is proposed to approximate second-order spatial derivatives directly in this paper to eliminate the unphysical oscillations.In this method, the second-order Taylor series is used, and symmetric moment matrices are obtained.It can directly approximate second-order derivatives.Numerical characteristics including truncation error, numerical accuracy, convergence rate, and stability of both FO-SSPH and SO-SSPH are analyzed theoretically and numerically.One key superiority of SO-SSPH to the MSPH method [17] and the scheme of Fatehi and Manzari [30] is that no kernel derivatives are needed in the calculation and the moment matrix is symmetric.This offers many numerical advantages including more choices for the kernel and saving computing resources.Numerical results show that, for transient heat conduction with initial smooth distribution, both FO-SSPH and SO-SSPH can achieve second order convergence rate, which is consistent with the theoretical analysis.However, for both one-and two-dimensional transient heat conduction problems with a discontinuous initial distribution, the FO-SSPH method has only a first-order convergence rate due to numerical oscillations, while the SO-SSPH method avoids unphysical oscillations and achieves second-order convergence.

Transient Heat Conduction Problem
In the Cartesian coordinate system, the governing equation for two-dimensional transient heat conduction problems reads where T is temperature, ρ is material density, c p is material specific heat capacity, G is source term, Ω is the computational domain, k x and k y are thermal conductivity coefficients in x and y directions, respectively.These coefficients are influenced by many factors, such as temperature and material.
For simplicity, we assume that both of them are constants, meaning that the transient heat conduction equation is linear.
To close Equation (1), the initial condition is and Dirichlet boundary condition as follows and Neumann boundary condition as follows where T 0 is initial temperature distribution, T 1 is boundary temperature, bT is heat source on the boundary, q = (q x , q y ) = (−k x ∂T ∂x , −k y ∂T ∂y ) represents heat flux, n = (n x , n y ) T is the unit normal vector on the boundary, Γ D ∪ Γ N = ∂Ω is the boundary of the domain, and If we further assume that there is no heat source, then Equation (1) reduces to where κ x = k x /ρc p and κ y = k y /ρc p .Introducing the heat fluxes as intermediate variables, the second-order heat conduction Equation ( 5) can be transformed into the following first-order system as

SPH Method
In SPH, the domain Ω is represented by a finite number of particles, each of which has density, mass, temperature, and other relevant physical quantities.Kernel approximation and particle approximation are the two key ideas in standard SPH method [6,7].
For any physical quantity f at position x i = (x i , y i ), the SPH kernel approximation of its gradient can be expressed as [7] ∇ where ∇W is the kernel gradient, and h is the smoothing length.The kernel is required to own properties of normalization, symmetry, compactness, and convergence to Dirac function [7].The cubic spline function is usually employed as the kernel.
Applying particle approximation to the integral terms in Equation ( 7), we obtain where and ∆V j is the volume of jth particle.

SO-SSPH Method
The Taylor series for a function f (x, y) around point (x i , y i ) in a two-dimensional domain gives where We neglect the third-and higher-order derivatives and multiply both sides with a compactly supported function Then, integrating over the computational domain and approximating the integrals with Riemman sum, we can get a system with five unknowns in matrix form as where Assume that W k = P k W ji and P 1 = x ji , P 2 = y ji , P 3 = 1 2 x ji 2 , P 4 = 1 2 y ji 2 , P 5 = x ji y ji , where x ji = x j − x i , y ji = y j − y i .Then the second derivatives are approximated as As the function is expanded up to second derivatives, and the moment matrix M 5×5 is symmetric, it is referred to as second-order symmetric smoothed particle hydrodynamics (SO-SSPH) method.
When the function is expanded up to first derivatives, the resulting system becomes where The first derivatives are approximated as This was referred to as the first-order symmetric smoothed particle hydrodynamics (FO-SSPH) method by Jiang et al. [32].
To effectively solve heat conduction problems with discontinuous initial distributions, different strategies for spatial derivative approximations are investigated in this paper.Approximating the second-order spatial derivatives in Equation ( 5) using the SO-SSPH method, and approximating the heat flux terms and their derivatives in Equation ( 6) using the FO-SSPH method were carried out.In theory, time integration can be completed by any effective algorithm, and for convenience, we use the Euler forward difference scheme herein.

Truncation Error
To approximate the spatial derivatives, both FO-SSPH and SO-SSPH methods reserve the first several terms and neglect the rest in Taylor series expansion.Obviously, this introduces deviation, which is called truncation error [34,35] and denoted by TE.
In the SO-SSPH method, the exact matrix form for the derivative approximations gives the truncation error corresponding to the approximation of second-order spatial derivatives in SO-SSPH method writes where β x and β y are the coefficients of ∂y 2 , respectively.In the FO-SSPH method, the exact matrix form for the derivatives gives Therefore, the truncation error corresponding to the approximation of first-order derivatives writes where α x and α y are the coefficients of

Numerical Accuracy
According to the definition of matrix M 5×5 in Equation (10), every matrix element can be represented as ∑ j x α ji y β ji W ji ∆V j , (α, β ∈ N + ).Considering the symmetry property of the kernel function, the following conclusion can be drawn for evenly distributed particles as applied in this work.
As the magnitudes of particle distances are the same, we obtain where h = (∆x, ∆y) T .Therefore, the moment matrix satisfies It can be seen that M 5×5 belongs to a partitioned diagonal matrix, so its inverse matrix follows as Similarly, Therefore, Equation ( 15) can be expressed as That is to say, using the SO-SSPH method to approximate second derivatives possesses second-order accuracy.
According to the definition of matrix M 2×2 in Equation ( 12) and similar analysis, we obtain Therefore, Equation ( 17) can be expressed as That is to say, using the FO-SSPH method to approximate first derivatives possesses second-order accuracy too.

Convergence Rate
The convergence is theoretically defined to describe the degree of the closeness between numerical solution and exact solution when (∆t, h) → (0, 0) [36].The actual convergence rates not only depend on the discretization procedure, but also on the smoothness of the solution [37].
Assuming that the numerical error satisfies and after taking the logarithm, we obtain log E = q log|h|+ log C Therefore, the function in Equation ( 28) represents a straight line in log-log coordinate system, and its slope describes the convergence rate of the corresponding numerical method.

Stability Analysis
Equation (11) gives the SO-SSPH discrete form to approximate second derivatives, which can also be written as with Substituting Equation ( 29) for spatial derivatives and forward difference for time derivative in Equation ( 5), we obtain where i is the imaginary unit and ν is the column vector of the wave number in the Fourier analysis method [34,35], h j = x j − x i T = ∆x j , ∆y j T .Thus, the amplification factor is Taking Equation ( 30) into consideration, it can be simplified to According to the von Neumann stability criterion [38], if the discrete form is stable, the amplification factor must satisfy the condition of |G| ≤ 1.For simplicity, we choose its necessary condition |Re(G)| ≤ 1 in this paper, which leads to Equation ( 13) gives the FO-SSPH discrete form to approximate first derivatives, which also can be written as with Substituting Equation (35) for spatial derivatives and forward difference for time derivative in Equation ( 6), we obtain x,i + ∑ j λ j q n+1 x,j ) − (γ i q n+1 y,i + ∑ j γ j q n+1 y,j ) which can be simplified to According to the Fourier analysis, together with Equation ( 36), the amplification factor follows and further simplification leads to According to the von Neumann stability criterion [38], if the discrete form is stable, the amplification factor must satisfy the condition of |G| ≤ 1.For simplicity, its necessary condition |Re(G)| ≤ 1 is employed and we obtain

One-Dimensional Case
Considering the one-dimensional transient heat conduction problem with discontinuous initial condition and boundary condition The exact solution is in which er f is the error function defined as er f (x) = 2 √ π x 0 e −µ 2 dµ.The parameters of L = 2 m, a = 0.5 m, κ x = k x /ρc p = 0.01 m 2 /s, and T 0 = 10 • C are used in the simulation.The initial temperature distribution is given in Figure 1a.Both FO-SSPH and SO-SSPH use the same parameters, particle distance ∆x = 0.1 m, and smoothing length h = 1.05 ∆x.According to the stability analysis given above, the stability conditions of FO-SSPH and SO-SSPH are ∆t ≤ 2 s and ∆t ≤ 0.5 s, respectively.To avoid the influence of the time integration error on the subsequent spatial convergence analysis, a sufficiently small time step of ∆t = 0.00005 s is chosen, and the simulation time is 5 s.
The exact solution is  Numerical oscillations appearing in the FO-SSPH solution are mainly due to the discontinuity in the initial temperature.The temperature gradient at the discontinuous position tends to infinity in theory, while the approximated first derivative in Taylor expansion can only be a finite value.This Numerical results of FO-SSPH and SO-SSPH are shown in Figure 1b together with the exact solution.It is seen that FO-SSPH [32] suffers from unphysical oscillations, which do not exist in the SO-SSPH solution.Temperature distributions at four different times are shown in Figure 1c where excellent agreements with exact solutions are observed as expected.
Numerical oscillations appearing in the FO-SSPH solution are mainly due to the discontinuity in the initial temperature.The temperature gradient at the discontinuous position tends to infinity in theory, while the approximated first derivative in Taylor expansion can only be a finite value.This leads to degeneration of numerical accuracy.Therefore, for transient heat conduction problems with discontinuous initial distribution, the second-order spatial derivatives should be discretized directly to reduce error accumulation.
For further verification, the forward time central space (FTCS) method [33] is also used to solve this problem (first derivative is approximated by the central difference) and the results are shown as Figure 1d.It demonstrates that FTCS with heat fluxes suffers from unphysical oscillations too, which is consistent with the FO-SSPH results.However, without introducing the heat fluxes, FTCS (second derivative is approximated by the central difference) gives solutions without oscillations (see Figure 1d).
To further analyze the properties of the methods, we define the relative numerical error as where T e is exact temperature, T n is numerical temperature, and the summation is over all particles.
Varying the values of ∆x, six pairs of error data corresponding to FO-SSPH and SO-SSPH are calculated and the results are shown in Table 1.According to the data in Table 1, least-squares-based curve fitting is performed.Convergence rates in double logarithmic coordinates are shown in Figure 2a.It is found that FO-SSPH can only achieve a first-order convergence rate, while SO-SSPH has a second-order convergence rate.In Section 3, we have proven that both FO-SSPH for first derivative approximation and SO-SSPH for second derivative approximation possess second-order accuracy.It implies that the actual convergence rate depends not only on the accuracy of the method, but also on the smoothness of the solution.This is consistent with the Kuzmin's conclusion for mesh-based methods [37].leads to degeneration of numerical accuracy.Therefore, for transient heat conduction problems with discontinuous initial distribution, the second-order spatial derivatives should be discretized directly to reduce error accumulation.
For further verification, the forward time central space (FTCS) method [33] is also used to solve this problem (first derivative is approximated by the central difference) and the results are shown as Figure 1d.It demonstrates that FTCS with heat fluxes suffers from unphysical oscillations too, which is consistent with the FO-SSPH results.However, without introducing the heat fluxes, FTCS (second derivative is approximated by the central difference) gives solutions without oscillations (see Figure 1d).
To further analyze the properties of the methods, we define the relative numerical error as where e T is exact temperature, n T is numerical temperature, and the summation is over all particles.Varying the values of x Δ , six pairs of error data corresponding to FO-SSPH and SO-SSPH are calculated and the results are shown in Table 1.According to the data in Table 1, least-squares-based curve fitting is performed.Convergence rates in double logarithmic coordinates are shown in Figure 2a.It is found that FO-SSPH can only achieve a first-order convergence rate, while SO-SSPH has a second-order convergence rate.In Section 3, we have proven that both FO-SSPH for first derivative approximation and SO-SSPH for second derivative approximation possess second-order accuracy.It implies that the actual convergence rate depends not only on the accuracy of the method, but also on the smoothness of the solution.This is consistent with the Kuzmin's conclusion for mesh-based methods [37].For further verification, and to more clearly explain the FO-SSPH defects when dealing with heat conduction problems with discontinuous initial distribution, the performances of FO-SSPH and For further verification, and to more clearly explain the FO-SSPH defects when dealing with heat conduction problems with discontinuous initial distribution, the performances of FO-SSPH and SO-SSPH methods to solve a heat conduction problem with smooth (Gaussian) initial distribution are examined.The Gaussian initial distribution follows where M = √ 2πσ 0 and σ 0 = 2 m.The exact solution is written as The corresponding errors are shown in Table 2, and the convergence rate curves of FO-SSPH and SO-SSPH are drawn in Figure 2b.It is seen that FO-SSPH method can achieve second-order convergence, which is in accordance with the theoretical analysis.Therefore, when the FO-SSPH method is used to deal with heat conduction problems with discontinuous initial distributions, the theoretical convergence rate cannot be achieved due to numerical oscillations.

Two-Dimensional Case
Consider the two-dimensional transient heat Equation (5) with initial condition and boundary condition where a < L, b < L. The initial distribution is discontinuous and the exact solution is The FO-SSPH and SO-SSPH are used to simulate this problem and the results at cross-section y = 0.05 m are shown in Figure 3b together with the exact solution.It is seen that for two-dimensional heat conduction with discontinuous initial distribution, the FO-SSPH method still suffers from unphysical oscillations as might be the expected, but the proposed SO-SSPH does not.
To show the two-dimensional heat conduction process, temperature contours at five different times are depicted in Figure 4a-e.It is seen that with increasing time, high temperature area expands and amplitude decreases, which conforms to the law of heat conduction.In addition, each contour line is smooth, which demonstrates that there is no unphysical oscillations in the SO-SSPH solution as expected.This is more clear as shown in Figure 4f for the temperature distribution along y = 0.05 m.FO-SSPH and SO-SSPH are used to simulate this problem and the results at cross-section y = 0.05 m are shown in Figure 3b together with the exact solution.It is seen that for two-dimensional heat conduction with discontinuous initial distribution, the FO-SSPH method still suffers from unphysical oscillations as might be the expected, but the proposed SO-SSPH does not.
To show the two-dimensional heat conduction process, temperature contours at five different times are depicted in Figure 4a-e.It is seen that with increasing time, high temperature area expands and amplitude decreases, which conforms to the law of heat conduction.In addition, each contour line is smooth, which demonstrates that there is no unphysical oscillations in the SO-SSPH solution as expected.This is more clear as shown in Figure 4f for the temperature distribution along y = 0.05 m.
Processes 2018, 6, x FOR PEER REVIEW 13 of 18 FO-SSPH and SO-SSPH are used to simulate this problem and the results at cross-section y = 0.05 m are shown in Figure 3b together with the exact solution.It is seen that for two-dimensional heat conduction with discontinuous initial distribution, the FO-SSPH method still suffers from unphysical oscillations as might be the expected, but the proposed SO-SSPH does not.
To show the two-dimensional heat conduction process, temperature contours at five different times are depicted in Figure 4a-e.It is seen that with increasing time, high temperature area expands and amplitude decreases, which conforms to the law of heat conduction.In addition, each contour line is smooth, which demonstrates that there is no unphysical oscillations in the SO-SSPH solution as expected.This is more clear as shown in Figure 4f for the temperature distribution along y = 0.05 m.Varying the values of ∆x and ∆y, four pairs of error data corresponding to FO-SSPH and SO-SSPH are listed in Table 3.The convergence rate curves of FO-SSPH and SO-SSPH are shown in Figure 5.It is clear that similar to the one-dimensional case, FO-SSPH only achieves first-order convergence, whereas SO-SSPH has second-order convergence.Varying the values of ∆x and ∆y, four pairs of error data corresponding to FO-SSPH and SO-SSPH are listed in Table 3.The convergence rate curves of FO-SSPH and SO-SSPH are shown in Figure 5.It is clear that similar to the one-dimensional case, FO-SSPH only achieves first-order convergence, whereas SO-SSPH has second-order convergence.Varying the values of ∆x and ∆y, four pairs of error data corresponding to FO-SSPH and SO-SSPH are listed in Table 3.The convergence rate curves of FO-SSPH and SO-SSPH are shown in Figure 5.It is clear that similar to the one-dimensional case, FO-SSPH only achieves first-order convergence, whereas SO-SSPH has second-order convergence.

Discussion and Conclusions
For transient heat conduction problems with discontinuous initial distribution, a second-order symmetric smoothed particle hydrodynamics (SO-SSPH) method is proposed to approximate the second-order spatial derivative terms directly.Compared with the simulation results of the FO-SSPH method, the SO-SSPH method can effectively eliminate the unphysical oscillations in FO-SSPH.This is mainly due to the initial temperature which is discontinuous such that the local gradient tends to infinity, while the approximated first derivative can only be of a finite value.This implies that FO-SSPH is incapable of approximating the second derivatives indirectly in essence.In addition, according to the convergence analysis, it is found that for transient heat conduction problems with smooth initial distributions, both FO-SSPH and SO-SSPH achieve second-order convergence rates.However, for heat conduction problems with discontinuous initial distributions, FO-SSPH can only achieve first-order convergence rate (rather than the theoretical second-order accuracy) due to numerical oscillations, while our SO-SSPH method can achieve a second-order convergence rate.The actual convergence rates of the meshless methods FO-SSPH and SO-SSPH depend not only on the accuracy of the method, but also on the smoothness of the solution.The proposed SO-SSPH method can be applied to heat conduction problems with and without discontinuous initial distributions.Acknowledgments: We thank Arris Tijsseling at Eindhoven University of Technology, The Netherlands, for the inspiring discussion and useful suggestions.The authors would also like to express their sincere gratitude to the anonymous reviewers and the editor for many valuable suggestions and comments.

Conflicts of Interest:
The authors declare no conflicts of interest.The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

) in which erf is the error function defined as 1 .
T = ℃ are used in the simulation.The initial temperature distribution is given in Figure1a.Both FO-SSPH and SO-SSPH use the same parameters, particle distance According to the stability analysis given above, the stability conditions of FO-SSPH and SO-SSPH are avoid the influence of the time integration error on the subsequent spatial convergence analysis, a sufficiently small time step of and the simulation time is 5 s.Numerical results of FO-SSPH and SO-SSPH are shown in Figure1btogether with the exact solution.It is seen that FO-SSPH[32] suffers from unphysical oscillations, which do not exist in the SO-SSPH solution.Temperature distributions at four different times are shown in Figure1cwhere excellent agreements with exact solutions are observed as expected.

Figure 1 .
Figure 1.One-dimensional heat conduction with discontinuous initial temperature distribution: (a) initial condition; (b) results of FO-SSPH and SO-SSPH compared with exact solution; (c) SO-SSPH solutions at difference times; and (d) FTCS solutions for Equation (6) (with heat fluxes) and Equation (5) (without heat fluxes).

Figure 1 .
Figure 1.One-dimensional heat conduction with discontinuous initial temperature distribution: (a) initial condition; (b) results of FO-SSPH and SO-SSPH compared with exact solution; (c) SO-SSPH solutions at difference times; and (d) FTCS solutions for Equation (6) (with heat fluxes) and Equation (5) (without heat fluxes).

Figure 2 .
Figure 2. Convergence rate of FO-SSPH and SO-SSPH for solving one-dimensional heat conduction with (a) discontinuous initial distribution and (b) smooth initial distribution (Gaussian).

Figure 2 .
Figure 2. Convergence rate of FO-SSPH and SO-SSPH for solving one-dimensional heat conduction with (a) discontinuous initial distribution and (b) smooth initial distribution (Gaussian).

Figure 3 .
Figure 3. Two-dimensional heat conduction with discontinuous temperature distribution: (a) initial condition and (b) results of FO-SSPH and SO-SSPH along y = 0.05 m compared with exact solution.

Figure 3 .
Figure 3. Two-dimensional heat conduction with discontinuous temperature distribution: (a) initial condition and (b) results of FO-SSPH and SO-SSPH along y = 0.05 m compared with exact solution.

Figure 3 .
Figure 3. Two-dimensional heat conduction with discontinuous temperature distribution: (a) initial condition and (b) results of FO-SSPH and SO-SSPH along y = 0.05 m compared with exact solution.

Figure 5 .
Figure 5. Convergence rate of FO-SSPH and SO-SSPH for solving two-dimensional heat conduction with discontinuous initial distribution.

Figure 5 .
Figure 5. Convergence rate of FO-SSPH and SO-SSPH for solving two-dimensional heat conduction with discontinuous initial distribution.

Figure 5 .
Figure 5. Convergence rate of FO-SSPH and SO-SSPH for solving two-dimensional heat conduction with discontinuous initial distribution.

Funding:
This research was funded by the National Natural Science Foundation of China (Grant No. 51478305) and the National Key Research and Development Program of China (Grant No. 2016YFC0401900, Grant No. 2018YFC1407400).

Table 1 .
Errors of FO-SSPH and SO-SSPH for one-dimensional heat conduction with discontinuous initial distribution.

Table 1 .
Errors of FO-SSPH and SO-SSPH for one-dimensional heat conduction with discontinuous initial distribution.

Table 3 .
Errors of FO-SSPH and SO-SSPH for two-dimensional heat conduction with discontinuous initial distribution.

Table 3 .
Errors of FO-SSPH and SO-SSPH for two-dimensional heat conduction with discontinuous initial distribution.

Table 3 .
Errors of FO-SSPH and SO-SSPH for two-dimensional heat conduction with discontinuous initial distribution.