# Study on SPH Viscosity Term Formulations
## Materials and Methods I--General ISPH Methodology  
The basic concept of SPH is that a particle is fundamental in the Lagrangian description and the motion of a continuum can be represented with high accuracy by simulating the motion of a large number of such particles. Through the use of the integral interpolants, any field variables can be expressed by the integrals that are approximated by the summation interpolants over the neighboring particles. Similarly, the spatial derivatives at each particle location can also be evaluated by the summation interpolants. The derivatives are computed based on the kernel function, which can be obtained analytically. For example, the particle approximation of a value \(f(\mathbf{r}_{i})\) and its gradient \(\
abla f(\mathbf{r}_{i})\) can be shown as  
\[f(\mathbf{r}_{i})=\sum_{j=1}^{N}\frac{m_{j}}{\rho_{j}}f(\mathbf{r}_{j})W( \mathbf{r}_{i}-\mathbf{r}_{j},h_{i}) \tag{1}\]  
\[\
abla f(\mathbf{r}_{i})=\sum_{j=1}^{N}\frac{m_{j}}{\rho_{j}}f(\mathbf{r}_{j}) \
abla_{i}W(\mathbf{r}_{i}-\mathbf{r}_{j},h_{i}) \tag{2}\]  
where \(i\) is the reference particle and \(j\) represents all the neighboring particles; \(N\) is the number of particles for the summation; \(m_{j}\) is the particle mass; \(\rho_{j}\) is the particle density; \(\mathbf{r}_{i}\) and \(\mathbf{r}_{j}\) are the positions of the particle; \(W\) is the SPH interpolation kernel function; \(\
abla_{i}W\) represents the gradient of the kernel taken with respect to the reference particle; and \(h\) is the kernel smoothing length. The kernel function can have many different forms and a cubic spline kernel is used in this paper due to its computational accuracy and efficiency. In a practical SPH computation, the gradient of the pressure \(P\) and divergence of the velocity \(\mathbf{u}\) are usually calculated by  
\[\
abla P_{i}=\rho_{i}\sum_{j=1}^{N}m_{j}\left(\frac{P_{j}}{\rho_{j}^{2}}+\frac {P_{i}}{\rho_{i}^{2}}\right)\
abla_{i}W(\mathbf{r}_{i}-\mathbf{r}_{j},h_{i}) \tag{3}\]  
\[\
abla\cdot\mathbf{u}_{i}=-\frac{1}{\rho_{i}}\sum_{i=1}^{N}m_{j}\mathbf{u}_{ij }\cdot\
abla_{i}W(\mathbf{r}_{i}-\mathbf{r}_{j},h_{i}) \tag{4}\]  
where \(\mathbf{u}_{ij}=\mathbf{u}_{i}-\mathbf{u}_{j}\) is defined as the velocity difference. In an ISPH approach the fluid density is assumed to be constant. Therefore, the mass and momentum conservation equations, i.e., N-S equations, can be written in the following quasi-Lagrangian description as  
\[\
abla\cdot\mathbf{u}=0 \tag{5}\]  
\[\frac{D\mathbf{u}}{Dt}=-\frac{1}{\rho}\
abla P+\mathbf{g}+v_{0}\
abla^{2} \mathbf{u} \tag{6}\]  
where \(t\) is the time; \(\mathbf{g}\) is the gravitational acceleration; and \(v_{0}\) is the kinematic viscosity. The two-step projection method [24] is used to solve the velocity and pressure in Equations (5) and (6). The first step is the prediction of velocity in the time domain without including the pressure term. The intermediate particle velocity and position are obtained as  
\[\mathbf{u}_{*}=\mathbf{u}_{t}+\Delta\mathbf{u}_{*} \tag{7}\]  
\[\Delta\mathbf{u}_{*}=\left(\mathbf{g}+v_{0}\
abla^{2}\mathbf{u}\right)\Delta t \tag{8}\]  
\[\mathbf{r}_{*}=\mathbf{r}_{t}+\mathbf{u}_{*}\Delta t \tag{9}\]  
where \(\mathbf{u}_{t}\) and \(\mathbf{r}_{t}\) are the velocity and position of particle at time \(t\); \(\Delta t\) is the time step; \(\Delta\mathbf{u}_{*}\) is the velocity change; and \(\mathbf{u}_{*}\) and \(\mathbf{r}_{*}\) are the intermediate velocity and position of particle.