## 2 Simulation Methodology
### Basic Theory of the Smoothed Particle Hydrodynamics Method
#### 2.1.2 Integral Interpolation Theory  
The SPH particles interact with each other by a kernel function, so two steps are crucial to the construction of the SPH equation, that is, kernel function interpolation and particle approximation. The function at a point in the domain \(\Omega\) is defined as the field function, which is expressed as follows:  
\[f(r)=\int f(r^{\prime})W(r-r^{\prime},h)dr^{\prime} \tag{7}\]  
where \(r\) is the spatial position vector; \(W(r-r^{\prime},h)\) is the kernel function, which depends on the distance between two particles \(|r-r^{\prime}|\) and the smoothing length; and the smoothing factor \(\kappa\) and \(h\) decide the affected domain of the kernel function (as shown in Figure 1). The spatial derivative of the field function is  
\[\
abla\cdot f(r)=-\int f(r^{\prime})\cdot\
abla W(r-r^{\prime},h)dr^{\prime} \tag{8}\]  
By virtue of integration, the kernel function interpolation delivers the field function and its spatial derivative, and then the particle approximation is conducted to discretize the kernel function interpolation. In a certain domain of solution, the integral of the kernel function interpolation can be transformed through a set of particles into the discrete form. The particle mass \(m_{j}\) can be expressed as  
\[m_{j}=\rho_{j}\Delta V_{j}\left(j=1,2,\cdots,N\right) \tag{9}\]  
Figure 1: Affected domain of the kernel function. where \(m_{j}\)\(\rho_{j}\) and \(\Delta V_{j}\) are the mass, density, and volume of the particle \(j\), respectively; \(N\) is the total number of particles in the affected domain. After discretizing Equations (7) and (8), we obtain  
\[f(r_{i})=\sum_{j=1}^{N}\frac{m_{j}}{\rho_{j}}f(r_{j})W_{ij} \tag{10}\]  
\[\
abla\cdot f(r_{i})=\sum_{j=1}^{N}\frac{m_{j}}{\rho_{j}}f(r_{j})\
abla_{i}W_{ ij} \tag{11}\]  
\[\
abla_{i}W_{ij}=\frac{r_{ij}}{\left|r_{i}-r_{j}\right|}\frac{\partial W_{ij}} {\partial\left|r_{i}-r_{j}\right|} \tag{12}\]  
where \(W_{ij}\) is the abbreviation for kernel function \(W(r_{i}-r_{j},h)\). Equation (11) indicates that the field function value at the position of particle \(i\) can be obtained by weighting the kernel function values of all particles in the domain. Equation (12) reveals that the spatial derivative value of the field function depends on the kernel function gradient and the function values of particles in the affected domain. With Equations (10) and (11), the integral form of a function and its spatial derivative is transformed into the discrete summation of particles so that the SPH method can achieve numerical computation without mesh. In order to improve interpolation accuracy and computational stability, this paper uses the Moving Least Squares (MLS) method to correct the integration scheme [35; 36]. The kernel function chosen in this paper is the cubic spline kernel function. \[W_{ij}=A^{\prime}\times \begin{cases}1-\frac{3}{2}q^{2}+\frac{3}{4}q^{3}&0\leq q <1\\ \frac{1}{4}(2-q)^{3}&1\leq q<\kappa\\ 0&q\geq\kappa\end{cases}  \tag{13}\]  
where \(q=\left|r_{i}-r_{j}\right|/h\), \(A^{\prime}\) is taken \(\frac{1}{\pi h^{3}}\). The traditional SPH method suffers from tensile instability; this paper adopts the artificial stress method to solve this problem [37]. \[\frac{dv_{i}}{dt}=-\sum_{j=1}^{N}m_{j}(\frac{P_{i}+P_{j}}{\rho_{i}\rho_{j}}+f_ {ij}^{\psi}R_{ij})\
abla_{i}W_{ij}+\sum_{j}^{N}m_{j}\frac{\eta_{i}+\eta_{j}}{ \rho_{i}\rho_{j}}+\sum_{j}^{N}m_{j}\frac{\eta_{i}+\eta_{j}}{\rho_{i}\rho_{j}} v_{ij}(\frac{1}{r_{ij}}\frac{\partial W_{ij}}{\partial r_{ij}})+g \tag{14}\]  
where \(P_{i}\), \(P_{j}\) are the pressures of particle \(i\), \(j\). \(\eta_{i}\), \(\eta_{j}\) are dynamic viscosity of particle \(i\), \(j\). \(f_{ij}=W_{ij}/W(\Delta p,h)\), \(\Delta p\) is the initial distance between particles. \(\psi\) is the variable exponent based on the kernel function. \(v_{ij}=v_{i}-v_{j}\) denotes relative velocity. \(R_{ij}\) is an artificial stress term. \[R_{ij}=R_{i}+R_{j} \tag{15}\]  
\(R_{i}\) and \(R_{j}\) are artificial parameters determined by pressure and density. \[