## 3 SPH Formulations
### SPH Implementation Algorithm for Incompressible Fluid  
There are two ways to ensure the incompressibility of fluid in SPH simulations. One is to divide the momentum equation into two parts. In the first part, the momentum equation is solved in the presence of body forces and viscosity forces. The second-part equation is solved in the presence of body force and viscosity forces. The second part considers the effect of the pressure and ensures incompressibility for the mediums. The other method isemployed in this paper, including two steps as well [10]. In the former ISPH method, an intermediate density needs to be calculated according to the intermediate velocity values. However, in the present incompressible SPH method, constant density is used in the whole process hence the computation time is saved and no consideration about density oscillation is required. The whole process of the algorithm is shown in Section 5.  
\[\mathbf{v}^{n+1/2}=\mathbf{v}^{n}+\mathbf{g}\Delta t+\left(\frac{1}{\rho} \
abla\cdot\mathbf{v}^{n}\right)\Delta t \tag{15}\]  
For Newtonian fluids such as water, the viscosity coefficient \(\mu_{eff}\) (effective viscosity) has a constant value of \(\mu\). Hence, the SPH formulation of viscosity term can be written as [8]:  
\[\left(\frac{1}{\rho}\
abla\cdot\mathbf{\tau}\right)_{a}=\left(\sum\limits_{b} \frac{4m_{b}\mu}{(\rho_{a}+\rho_{b})^{2}}\frac{\mathbf{r}_{ab}\cdot(\mathbf{v} _{a}-\mathbf{v}_{b})\
abla W}{\left|\mathbf{r}_{ab}\right|^{2}+\eta^{2}}\right) \tag{16}\]  
where \(\mathbf{r}_{ab}\) means \(\mathbf{r}_{a}-\mathbf{r}_{b}\).  
\[\frac{\mathbf{v}^{n+1}-\mathbf{v}^{n+1/2}}{\Delta t}=-\frac{1}{\rho}\
abla P^ {n+1} \tag{17}\]  
Taking the divergence of Equation (17) and substituting it into Equation (5),  
\[\Delta P^{n+1}=\frac{\rho}{\Delta t}\
abla\cdot\mathbf{v}^{n+1/2} \tag{18}\]  
Equation (18) can be expressed into SPH form:  
\[\sum\limits_{b=1}^{N}\frac{m_{b}p_{ab}^{n+1}\mathbf{r}_{ab}^{n}}{\overline{ \mathbf{r}_{ab}^{2}+\eta^{2}}}\
abla W_{ab}^{m}=-\frac{\rho}{2\Delta t}\sum \limits_{b=1}^{N}m_{b}p_{ab}^{n+1/2}\
abla W_{ab}^{n} \tag{19}\]  
where \(n\) indicates a time step; \(\eta\) is a random small number to avoid zero denominator. The pressure \(P\) of particle \"\(a\)\" can be updated by Equation (19) which can be solved efficiently by the BI-CGSTAB method. Solid particles which are the neighboring particles of fluid particles are involved in the [PERSON]'s equation and pressure values are calculated simultaneously with fluid particles.  
Finally, the velocity and position of each fluid particle can be renewed for the next time step:  
\[\begin{array}{ll}\mathbf{v}^{n+1}&=\mathbf{v}^{n+1/2}-\frac{1}{\rho}\
abla P ^{n+1}\Delta t\\ &=\mathbf{v}^{n+1/2}-\sum\limits_{b=1}^{N}m_{b}\left(\frac{P^{n+1}}{\rho_{a}^ {2}}+\frac{P^{n+1}}{\rho_{b}^{2}}+\Pi_{ab}\right)\cdot\
abla W(\mathbf{x}_{a} -\mathbf{x}_{b},h)\end{array} \tag{20}\]  
where \(\Pi_{ab}\) is artificial viscosity. It is added into the momentum equation to avoid unphysical particle oscillation [5]:  
\[\prod\
olimits_{ij}=\left\{\begin{array}{l}\frac{-\Delta t_{ab}\mu_{ab}+ \beta\mu_{ab}^{2}}{\overline{\mu}_{ab}};\dots\mathbf{v}_{ab}\cdot\mathbf{r}_{ ab}<0\\ 0;\dots\dots\dots\dots\dots\dots\dots\mathbf{v}_{ab}\cdot\mathbf{r}_{ab}>0 \end{array}\right. \tag{21}\]  
where \(\mu_{ab}=\frac{|\mathbf{r}_{ab}\cdot\mathbf{r}_{ab}}{\overline{\mathbf{r}_{ab }^{2}+\eta^{2}}}\). Any value \(\overline{f}_{ab}\) denotes \(\frac{f_{a}+f_{b}}{2}\).  
For fluid particles near the boundary, the velocity calculation is affected by the pressure values of the solid particles, so that the fluid particles could not penetrate the solid boundary:  
\[\mathbf{x}^{n+1}=\mathbf{x}^{n}+\mathbf{v}^{n+1}\Delta t \tag{22}\]