2. Kernel Approximation and Kernel Function Conditions in the SPH Method (Preparation for Spatial Discretization) | Growing the particle method, and its present state

This section explains a spatial discretization scheme in the Smoothed Particle Hydrodynamics (SPH) method. In the SPH method, a continuous physical quantity is approximated by a kernel function, and then the kernel approximation becomes a discretize form with particle quantities in the vicinity of the target point. For example, a scalar function \( \phi\) of the field defined on a three-dimensional (or two-dimensional) region \( \Omega\) is expressed in the form of a volume integral using the Dirac delta function \( \delta\) at an arbitrary point \( \boldsymbol{x}\) in the region.

$$ \phi(\boldsymbol{x})=\int_{\Omega} \phi(\boldsymbol{\xi}) \delta(\boldsymbol{\xi}-\boldsymbol{x}) \mathrm{d} \boldsymbol{\xi} $$

Here, the symbol \(\int_{\Omega}(\cdot) \mathrm{d} \boldsymbol{\xi}\), which denotes the volume integral, refers to the integration of the volume, where the integrant is a function of the variable \(\boldsymbol{\xi}\). In addition, the Dirac delta function satisfies

$$ \delta(\boldsymbol{x})=\left\{\begin{array}{ll} \infty, & \boldsymbol{x}=\mathbf{0} \\ 0, & \boldsymbol{x} \neq \mathbf{0} \end{array}\right. $$

and when a constant function \(\phi(\boldsymbol{x})=1\) is considered, this function has the following property:

$$ \int_{\Omega} \delta(\boldsymbol{\xi}) \mathrm{d} \boldsymbol{\xi}=1 $$

The \( \delta\) function can be regarded as a very cuspidal function,^*1 where the function value becomes \( \infty\) or zero at the origin. In the SPH method, instead of this cuspidal function, a smooth function is used for approximation.

$$ \phi(\boldsymbol{x})=\int_{\Omega} \phi(\boldsymbol{\xi}) \delta(\boldsymbol{\xi}-\boldsymbol{x}) \mathrm{d} \boldsymbol{\xi} $$ $$ \approx \int_{\Omega} \phi(\boldsymbol{\xi}) w_{h}^{\mathrm{SPH}}(|\boldsymbol{\xi}-\boldsymbol{x}|) \mathrm{d} \boldsymbol{\xi}=: \phi^{\mathrm{SPH}}(\boldsymbol{x}) $$

^*1 It is a unique function not only in its shape but also in mathematic properties, and is called a generalized function.

\(\phi^{\mathrm{SPH}}(\boldsymbol{x})\) is described as kernel approximation of \(\phi(\boldsymbol{x})\), and this smooth function \(w_{h}^{\mathrm{SPH}}\) is called the “kernel function.” It is assumed that \(h\) indicates a region (radius) to determine the target for smoothing, and in the SPH method, from the functions that satisfy the conditions below, a kernel function is selected.

$$ \text { (Compact support condition) } \quad w_{h}^{\mathrm{SPH}}(|\boldsymbol{\xi}-\boldsymbol{x}|)=0, \quad|\boldsymbol{\xi}-\boldsymbol{x}| \geq h $$ $$ \text { (Unity condition) }\quad \int_{\Omega} w_{h}^{\mathrm{SPH}}(|\boldsymbol{\xi}-\boldsymbol{x}|) d \boldsymbol{\xi}=1 $$

The above approximation Eq. (5) approximates the integral representation Eq. (4) using the delta function, and the compact support condition (6) and the unity condition (7) correspond to property (2) and (3) of the Dirac delta function, respectively. In the original integral representation, only the function value of the observation point \( \boldsymbol{x}\) is extracted from the dummy variable \(\boldsymbol{\xi}\), and therefore, the delta function \(\delta(\boldsymbol{\xi}-\boldsymbol{x}) \) has been utilized. In contrast, in the kernel approximation in the SPH method, points in the neighborhood of the observation point are also incorporated for smoothing, which means that the weighted average is obtained using the values in the neighborhood \(\boldsymbol{\xi}\), where \( |\boldsymbol{\xi}-\boldsymbol{x}|\) denoting the relative distance from the observation point \( \boldsymbol{x}\) is within a fixed radius \(h\). The compact support condition represents that the target for the weighted average is limited to the neighborhood within the fixed radius \(h\). The weight for such smoothing is the kernel function \(w_{h}^{\mathrm{SPH}}(|\boldsymbol{\xi}-\boldsymbol{x}|) \), where weight is given according to relative distance. The unity condition is a necessary condition to converge the kernel function to the delta function when the radius of influence \(h\) is converged to zero. In other words, this means that as the radius of influence is reduced, the weight of each point increases to satisfy the unity condition, and if the radius of influence is converged to the observation point as the limit, it should have infinite weight—in other words, it should be the Dirac delta function.

A kernel approximation in the integral form using this smooth function is the basis of the SPH method, and also became the origin of the term Smoothed Particle Hydrodynamics (SPH). Now, using the Taylor expansion, it is examined whether the equation, in which the Delta function has been changed to a kernel function, sufficiently holds as an approximation. Based on the compact support condition, the dummy variable \(\boldsymbol{\xi}\) appearing in the integrand could be considered as a point in the neighborhood of the observation point \( \boldsymbol{x}\). When Taylor expansion is performed around the observation point, where the function value in the neighborhood of this observation point \(\phi(\boldsymbol{\xi})=\phi(\boldsymbol{x}+\Delta \xi)\) is fixed,^*2 the following equations are obtained:

$$ \phi(\boldsymbol{\xi})=\phi(\boldsymbol{x})+\frac{\partial \phi}{\partial \boldsymbol{\xi}}(\boldsymbol{x}) \cdot \Delta \boldsymbol{\xi}+\frac{1}{2} \Delta \boldsymbol{\xi} \cdot \frac{\partial^{2} \phi}{\partial \boldsymbol{\xi}^{2}}(\boldsymbol{x}) \Delta \boldsymbol{\xi}+\mathcal{O}\left(\Delta \boldsymbol{\xi}^{3}\right) $$ $$ =\phi(\boldsymbol{x})+(\boldsymbol{\xi}-\boldsymbol{x}) \cdot \frac{\partial \phi}{\partial \boldsymbol{\xi}}(\boldsymbol{x})+\mathcal{O}\left(\Delta \boldsymbol{\xi}^{2}\right) $$

^*2 The same observation point is used at all points in the neighborhood to perform Taylor expansion.

where \(\Delta \boldsymbol{\xi}=\boldsymbol{\xi}-\boldsymbol{x}\) and hereafter, the approximation shall be conducted up to the first-order term of \(\Delta \boldsymbol{\xi}\).

When this Taylor expansion Eq. (9) is substituted for Eq. (5), the following equations are obtained:

$$ \begin{equation} \begin{aligned} \phi^{\mathrm{SPH}}(\boldsymbol{x}) &=\int_{\Omega} \phi(\boldsymbol{\xi}) w_{h}^{\mathrm{SPH}}(|\boldsymbol{\xi}-\boldsymbol{x}|) \mathrm{d} \boldsymbol{\xi} \\ &=\int_{\Omega}\left\{\phi(\boldsymbol{x})+(\boldsymbol{\xi}-\boldsymbol{x}) \cdot \frac{\partial \phi}{\partial \boldsymbol{\xi}}(\boldsymbol{x})\right\} w_{h}^{\mathrm{SPH}}(|\boldsymbol{\xi}-\boldsymbol{x}|) \mathrm{d} \boldsymbol{\xi}+R \end{aligned} \end{equation} $$ $$ =\phi(\boldsymbol{x}) \int_{\Omega} w_{h}^{\mathrm{SPH}}(|\boldsymbol{\xi}-\boldsymbol{x}|) \mathrm{d} \boldsymbol{\xi}+\frac{\partial \phi}{\partial \boldsymbol{\xi}}(\boldsymbol{x}) \cdot \int_{\Omega}(\boldsymbol{\xi}-\boldsymbol{x}) w_{h}^{\mathrm{SPH}}(|\boldsymbol{\xi}-\boldsymbol{x}|) \mathrm{d} \boldsymbol{\xi}+R $$ $$ =\phi(\boldsymbol{x}) \int_{\Omega} w_{h}^{\mathrm{SPH}}\left(\left|\boldsymbol{\xi}^{\prime}\right|\right) \mathrm{d} \boldsymbol{\xi}^{\prime}+\frac{\partial \phi}{\partial \boldsymbol{\xi}}(\boldsymbol{x}) \cdot \int_{\Omega} \boldsymbol{\xi}^{\prime} w_{h}^{\mathrm{SPH}}\left(\left|\boldsymbol{\xi}^{\prime}\right|\right) \mathrm{d} \boldsymbol{\xi}^{\prime}+R $$

\(R\) is an error related to \(\mathcal{O}\left(\Delta \boldsymbol{\xi}^{2}\right)\) in Eq. (9). In the deformation from Eq. (10) to Eq. (11), since the integrand is considered the function of the dummy variable \(\boldsymbol{\xi}\), \(\phi(\boldsymbol{x})\) and \(\frac{\partial \phi}{\partial \boldsymbol{\xi}}(\boldsymbol{x})\) are values at the fixed observation point \(\boldsymbol{x}\), and it is possible to move them outside the integration. Moreover, in Eq. (11) and Eq. (12), \(\boldsymbol{\xi}-\boldsymbol{x}\) denotes a relative position centering on the observation point \(\boldsymbol{x}\) and is expressed by \(\boldsymbol{\xi}^{\prime}\). Ultimately, to ensure the kernel approximation \(\phi^{\mathrm{SPH}}(\boldsymbol{x})\) shown in Eq. (11) holds as an approximation of the function \(\phi(\boldsymbol{x})\), the following two conditions need to be satisfied:

$$ \int_{\Omega} w_{h}^{\mathrm{SPH}}\left(\left|\boldsymbol{\xi}^{\prime}\right|\right) \mathrm{d} \boldsymbol{\xi}^{\prime}=1 $$ $$ \int_{\Omega} \boldsymbol{\xi}^{\prime} w_{h}^{\mathrm{SPH}}\left(\left|\boldsymbol{\xi}^{\prime}\right|\right) \mathrm{d} \boldsymbol{\xi}^{\prime}=0 $$

The first equation indicates that the value of the kernel function defined around the observation point \(\boldsymbol{x}\) should be “1” when it is integrated into the region; this is a condition for the kernel approximation referred to above as the “unity condition.” The second condition is satisfied if the kernel function is an even function (a function in which variables are symmetric with respect to the origin and function values are equal)

$$ w_{h}^{\mathrm{SPH}}\left(\left|\boldsymbol{\xi}^{\prime}\right|\right)=w_{h}^{\mathrm{SPH}}\left(\left|-\boldsymbol{\xi}^{\prime}\right|\right) $$

however, since the kernel function has been set as a function of distance \(r=|\boldsymbol{\xi}-\boldsymbol{x}|\), this property is automatically satisfied. Because of this kernel function property, \(\boldsymbol{\xi}^{\prime} w_{h}^{\mathrm{SPH}}\left(\left|\boldsymbol{\xi}^{\prime}\right|\right)\) becomes an odd function, indicating that the above condition Eq. (14) is satisfied. Lastly, to meet the demand to reduce the computation load as much as possible during a numerical simulation, i.e., by setting a narrow range for region integration, the compact support condition has been added. Now, all necessary conditions (and demands) given to the kernel function have been prepared. All the conditions are summarized below.

As the kernel function \(w_{h}^{\mathrm{SPH}}\) that satisfies the above seven conditions, piecewise polynomials called “cubic spline function” and “quintic spline function” are often used in the SPH numerical simulations in accordance with the policy to select the simplest possible functions (see Fig. 1).

$$ w_{h}^{\mathrm{SPHc}}(r)=\alpha_{d}^{\mathrm{c}}\left\{\begin{array}{ll} 1-6(r / h)^{2}+6(r / h)^{3}, & 0 \leq r / h<1 / 2, \\ 2(1-r / h)^{3}, & 1 / 2 \leq r / h<1, \\ 0, & r / h \geq 1, \end{array}\right. $$ $$ w_{h}^{\mathrm{SPHq}}(r)=\alpha_{d}^{\mathrm{q}}\left\{\begin{array}{ll} (1-r / h)^{5}-6(2 / 3-r / h)^{5}+15(1 / 3-r / h)^{5}, & 0 \leq r / h<1 / 3, \\ (1-r / h)^{5}-6(2 / 3-r / h)^{5}, & 1 / 3 \leq r / h<2 / 3, \\ (1-r / h)^{5}, & 2 / 3 \leq r / h<1, \\ 0, & r / h \geq 1 . \end{array}\right. $$

Here, \(\alpha_{d}^{c}\) and \(\alpha_{d}^{\mathrm{q}}\) are the constants to satisfy the unity condition (7), and can be determined in advance by setting the radius of influence as follows: \(\alpha_{d}^{c}=10 / 7 \pi h^{2} \) (two dimensions) or \(1 / \pi h^{3} \) (three dimensions); \(\alpha_{d}^{\mathrm{q}}=15309 / 478 \pi h^{2}\) (two dimensions) or \(2187 / 40 \pi h^{3} \) (three dimensions). Since the aforementioned spline functions are piecewise polynomials requiring case division, the Wendland function [8] requiring no case division is also sometimes used.