## 2 Numerical Method and Model Development  
### Mathematical Model Based on the SPH Method  
The SPH method is a Lagrangian particle method, wherein the computational domain is composed of a set of discrete particles. Each particle is assigned physical properties such as mass, density, pressure, and velocity. Based on these particles, the time derivatives of various physical quantities can be obtained using the SPH kernel approximation algorithm, which allows the time integration of the physical process of the flow. The kernel approximation of a field function \(f(\mathbf{r})\) at a spatial position \(\mathbf{r}\), \(\langle f(\mathbf{r})\rangle\), can be expressed as  
\[\langle f(\mathbf{r})\rangle=\int_{\Omega}f(\mathbf{r}^{*})W(\mathbf{r}-\mathbf{r}^{*},h)dV_{ \mathbf{r}}, \tag{1}\]  
where \(W(\mathbf{r}-\mathbf{r}^{*},h)\) is the smoothing kernel function [28], \(h\) is the smoothing length [29], and \(V\) is the volume. \(\Omega\) represents the support domain of the smoothing kernel function, which is a function of \(h\). A good tradeoff between the computational accuracy and the kernel approximation cost can be realized using the Wendland kernel function [30], which is expressed as \[W(\mathbf{r}-\mathbf{r}^{*},h)=\alpha_{D}\Big{(}1-\frac{q}{2}\Big{)}^{4}(2q+1) \tag{2}\]  
where \(D\) is the dimensions of the computational domain and \(q\) is the dimensionless variable, expressed as \(q=(r-r^{*})/h\). In addition, \(\alpha_{D}\) is a normalization factor which is given as \(7/(4\pi h^{2})\) for a two-dimensional (2D) domain and \(21/(16\pi h^{3})\) for a 3D domain. As the computational domain is composed of a set of discrete particles with particular physical properties, the continuous format of the kernel approximation, as expressed in Equation (1), should be transformed into a discrete one:  
\[\langle f(\mathbf{r}_{i})\rangle=\sum\limits_{j=1}^{N}f(\mathbf{r}_{j})W(\mathbf{r}_{i}- \mathbf{r}_{j},h)V_{j} \tag{3}\]  
where \(N\) is the total number of particles located in the support domain of particle \(i\). \(V_{j}\) is the virtual volume of the neighboring particle \(j\), which is the ratio of the mass \(m_{j}\) to the density \(\rho_{j}\). Consequently, the kernel approximation of the derivative of the field function in the discrete format, \(\langle\
abla\cdot f(\mathbf{r})\rangle\), can be derived as  
\[\langle\
abla\cdot f(\mathbf{r}_{i})\rangle=\sum\limits_{j=1}^{N}f(\mathbf{r}_{j}) \cdot\
abla_{i}W(\mathbf{r}_{i}-\mathbf{r}_{j},h)\frac{m_{j}}{\rho_{j}} \tag{4}\]  
Equation (4) indicates that the derivative of the field function can be transferred to the derivative of the smooth kernel function, thereby reducing the smoothness requirements for the field function and producing stable results in numerical approximations of PDEs [31]. Thus, the governing equations in the form of weakly compressible SPH are derived as  
\[\left\{\begin{aligned} \frac{d\rho_{i}}{dt}=-\rho_{i}\sum\limits_{j= 1}^{N}\left(\mathbf{u}_{i}-\mathbf{u}_{j}\right)\cdot\
abla_{i}W(\mathbf{r}_{i}-\mathbf{r}_{j},h)\frac{m_{j}}{\rho_{j}}\\ \frac{du_{i}}{dt}=\frac{1}{\rho_{i}}\sum\limits_{j=1}^{N}\left(P_ {{}_{i}}+P_{{}_{j}}\right)\
abla_{i}W(\mathbf{r}_{i}-\mathbf{r}_{j},h)\frac{m_{j}}{ \rho_{j}}+\mathbf{g}\\ \frac{d\mathbf{r}_{i}}{dt}=\mathbf{u}_{i}\\ P_{i}=\frac{\rho_{0}c_{0}^{2}}{\gamma}\bigg{(}\bigg{(}\frac{\rho_{i }}{\rho_{0}}\bigg{)}^{\gamma}-1\bigg{)}\end{aligned}\right. \tag{5}\]  
where \(\mathbf{u}\) is the velocity. The pressure of particle \(i\), \(P_{i}\), is expressed in the form of the Tait equation, where \(\rho_{0}\) is the density of the fluid when \(P=0\) Pa, \(c_{0}\) is the speed of sound [32; 33], and \(\gamma\) is a constant [34; 35] that usually takes the value of 7 in hydrodynamic problems. The continuity equation employs an asymmetric SPH approximation scheme to minimize errors arising from discontinuities [36], and the momentum equation adopts symmetric SPH approximation schemes to ensure better conservation [37].