## 3 Materials and Methods
### SPH Applied to Fluid Simulations  
In the Eulerian formulation, fluids are characterized by a velocity field \(\mathbf{v}\), a density field \(\rho\) and a pressure field \(p\). The progression of these fields over time is given by two equations. The first guarantees the conservation of mass  
\[\frac{\partial\rho}{\partial t}+\
abla\cdot(\rho\mathbf{v})=0, \tag{5}\]  
while the Navier-Stokes equation defines the conservation of momentum  
\[\rho(\frac{\partial v}{\partial t}+\mathbf{v}\cdot\
abla\mathbf{v})=-\
abla p +\rho\mathbf{g}+\mu\
abla^{2}\mathbf{v}, \tag{6}\]  
where \(\mathbf{g}\) is an external force and \(\mu\) the viscosity of the fluid.  
Using particles instead of a grid based formulation, Equation (6) can be expressed as  
\[\rho_{i}\mathbf{a}_{i}=\mathbf{f}_{i}^{\mathit{pressure}}+\mathbf{f}_{i}^{ \mathit{external}}+\mathbf{f}_{i}^{\mathit{viscosity}}, \tag{7}\]  
where \(\mathbf{a}_{i}\) corresponds to the acceleration of particle \(i\) and is integrated using the Leap-Frog scheme [54]:  
\[\mathbf{v}_{i+\frac{1}{2}}=\mathbf{v}_{i-\frac{1}{2}}+\Delta t \mathbf{a}_{i} \tag{8}\] \[\mathbf{r}_{i}=\mathbf{r}_{i-1}+\Delta t\mathbf{v}_{i-\frac{1}{2 }}.\]  
The time step size \(\Delta t\) that will be used must be adapted by the Courant-Freidrich-Levy condition [55]:  
\[\Delta t\leq 0.4\frac{d}{||v_{max}||}, \tag{9}\]  
where \(d\) is the particle diameter, and \(v_{max}\) is the maximum particle velocity.  
Substituting Equations (3) and (4) into the pressure and viscosity terms of the Navier-Stokes equation and solving according to [PERSON] et al. [53] yields:  
\[\mathbf{f}_{i}^{\mathit{pressure}}=-\sum_{j}m_{j}\frac{p_{i}+p_{j}}{2\rho_{j}} \
abla W(\mathbf{r}_{i}-\mathbf{r}_{j},h), \tag{10}\]  
\[\mathbf{f}_{i}^{\mathit{viscosity}}=\sum_{j}\frac{\mu_{i}+\mu_{j}}{2}m_{j}\frac {\mathbf{v}_{j}-\mathbf{v}_{i}}{\rho_{j}}\
abla^{2}W(\mathbf{r}_{i}-\mathbf{ r}_{j},h), \tag{11}\]  
where \(p\) is the pressure, \(\mathbf{v}\) the velocity, and \(\mu\) the viscosity coefficient. The pressure \(p_{i}\) of particle \(i\) is computed via the modified gas state equation suggested by [56]:  
\[p_{i}=k(\rho_{i}-\rho_{0}), \tag{12}\]  
where \(\rho_{0}\) is the rest density, and \(k\) is a stiffness constant that scales the pressure, and, thus, the pressure gradient and the respective pressure forces [34].  
SPH Smoothing Kernels for Fluid Simulations  
Stability, accuracy, and speed of the SPH method rely upon the choice of smoothing kernels. For fluid simulations [PERSON] et al. [53] designed the following kernel  
\[W_{\mathit{polyg}}(\mathbf{r},h)=\frac{315}{64\pi h^{9}}\begin{cases}(h^{2}-r^ {2})^{3}&0\leq r\leq h\\.0&\text{otherwise}\end{cases} \tag{13}\]  
Unfortunately, if this kernel is used for the calculation of the pressure forces, particles build clusters under high pressure, and as they get close to each other, the repulsion force vanishes because the gradient of the kernel approaches zero. [PERSON] [56] solved this issue by using a spiky kernel with a non vanishing gradient near the center. That kernel generates necessary repulsion forces in pressure calculations and is defined as:  
\[W_{\mathit{spiky}}(\mathbf{r},h)=\frac{15}{\pi h^{6}}\begin{cases}(h-r)^{3}& 0\leq r\leq h\\ 0&\text{otherwise}\end{cases}. \tag{14}\]  
Viscosity is caused by friction and, thus, should only have a smoothing effect on the velocity field. However, if a standard kernel is used for viscosity, the resulting viscosity forces do not always have this property. For this reason, for the computation of viscosity forces, [PERSON] et al. [53] proposed a third kernel:\[W_{viscosity}(\mathbf{r},h)=\frac{15}{2\pi h^{3}}\begin{cases}-\frac{r^{3}}{2h^{3}}+ \frac{r^{2}}{h^{2}}+\frac{h}{2r}-1&0\leq r\leq h\\ 0&\text{otherwise}\end{cases}. \tag{15}\]