To discretely solve the Navier-Stokes equations \[ \frac{\mathrm{d}\rho}{\mathrm{d}t} + \rho(\nabla \cdot \mathbf{v}) = 0, \quad \rho\left(\frac{\mathrm{D} \mathbf{v}}{\mathrm{D} t}\right) = - \nabla p + \mu \nabla^2 \mathbf{v} \] We can employ the Smoothed Particle Hydrodynamics method. This is a Lagrangian approach in fluid dynamics, where the fluid is represented by a finite set of particles. In this method, each particle indexed by \(i\) each carries a mass \(m_i\), position \(\mathbf{x}_i\), velocity \(\mathbf{v}_i\) and density \(\rho_i\). We use kernel interpolation to reconstruct quantities, \[ A(\mathbf{x}) = \sum_{j} m_j \frac{A_j}{\rho_j} W(\mathbf{x}-\mathbf{x}_j,h) \] where \(W\) is some radially symmetric smoothing kernel. The gradient of such quantity can also be discretized by \[ \nabla A(x) = \sum_{j} m_j \frac{A_j}{\rho_j}\nabla W(\mathbf{x}-\mathbf{x}_j,h) \] For our purposes, we use three kernels. For the calculation of most quantities, we use the Poly6 kernel. \[ W_{poly6}(\mathbf{r}, h) = \begin{cases} \dfrac{315}{64\pi h^9}(h^2 - \|\mathbf{r}\|^2)^3, & 0 \le \|\mathbf{r}\| \le h \\ 0, & \|\mathbf{r}\| > h \end{cases} \] The density \(\rho_i\) at particle \(i\) is then given by: \[ \rho_i = \sum_{j} m_j W_{poly6}(\mathbf{x}_i - \mathbf{x}_j, h) \] For pressure computations, we avoid the Poly6 kernel because its gradient vanishes at the origin. Instead, we employ the gradient of Debrun's spiky kernel: \[ \nabla W_{spiky}(\mathbf{r}, h) = \begin{cases} -\dfrac{45}{\pi h^6}(h - \|\mathbf{r}\|)^2 \dfrac{\mathbf{r}}{\|\mathbf{r}\|}, & 0 < \|\mathbf{r}\| \le h \\ 0, & \|\mathbf{r}\| > h \end{cases} \] Using the symmetric formulation to conserve momentum, the pressure force \(\mathbf{f}_i^p\) acting on particle \(i\) is: \[ \mathbf{f}_i^p = -\sum_{j} m_j \frac{p_i + p_j}{2\rho_j} \nabla W_{spiky}(\mathbf{x}_i - \mathbf{x}_j, h) \] where \(p_i = k(\rho_i - \rho_0)\) is the pressure derived from the ideal gas law. To compute the viscous shear forces, we require the Laplacian of the velocity field. A specialized viscosity kernel is used, where the Laplacian is defined as: \[ \nabla^2 W_{visc}(\mathbf{r}, h) = \begin{cases} \dfrac{45}{\pi h^6}(h - \|\mathbf{r}\|), & 0 \le \|\mathbf{r}\| \le h \\ 0, & \|\mathbf{r}\| > h \end{cases} \] The resulting viscosity force \(\mathbf{f}_i^v\) for particle \(i\) is calculated as: \[ \mathbf{f}_i^v = \mu \sum_{j} m_j \frac{\mathbf{v}_j - \mathbf{v}_i}{\rho_j} \nabla^2 W_{visc}(\mathbf{x}_i - \mathbf{x}_j, h) \] The final acceleration of a particle \(i\) is now calculated as \[ \mathbf{a}_i = \frac{\mathrm{D} \mathbf{v}_i}{\mathrm{D} t} = \frac{\mathbf{f}_i^p + \mathbf{f}_i^v}{\rho_i} + \mathbf{a}_{ext} \] Where \(\mathbf{a}_{ext}\) are external forces on the system. To numerically integrate the system, we use the leapfrog scheme. \[ \begin{aligned} \mathbf{v}_{i}^{t+1/2} &= \mathbf{v}_{i}^{t-1/2} + \mathbf{a}_i^{t} \Delta t\\ \mathbf{x}_{i}^{t+1} &= \mathbf{x}^t_i + \mathbf{v}_i^{t+1/2} \Delta t \end{aligned} \] Rendering was done by screen space fluid rendering.
Smoothed Particle Hydrodynamics
Screen Space Fluid Rendering