Assignment 3

Table of Contents

• 1 Getting Started
• 2 Implementation Details
  • 2.1 Density Computation
  • 2.2 Pressure Computation
  • 2.3 Acceleration Computation
  • 2.4 Boundary Conditions
  • 2.5 Artificial Viscosity
• 3 Interactive Controls
• 4 Debugging Tips
• 5 Submission
• 6 References

In this assignment, you will implement a 2D fluid simulation using the Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH) method, using the EUROGRAPHICS 2019 course notes by Koschier et al. as a reference for the details.

1 Getting Started

Pull the latest A3 starter code from: https://github.coecis.cornell.edu/cs5643/assignments-sp25. Running python3 main.py will display a GUI with initial particles.

The framework consists of four files:

• main.py: Simulation loop and GUI
• scene.py: Particle position initialization
• utils.py: SPH constants
• todo.py: Your implementation code goes here

To record a simulation video, use python3 main.py --video. Recording stops when you press the ESC key, generating both MP4 and GIF files. Here is the reference video. Note that this video is only for reference, and your results don’t need to be exactly the same.

2 Implementation Details

2.1 Density Computation

In SPH, a field $A$ is discretized using sampling points (particles) $j$ as:

$$⟨A(\overrightarrow{x})⟩=\sum _j{A}_j\frac{m_j}{{\rho }_j}W(\Vert \overrightarrow{x}-{\overrightarrow{x}}_j\Vert ,h)$$

Where $m$ is mass, $\rho$ is density, and $W$ is the kernel function. For density itself, the equation simplifies to:

$$⟨\rho (\overrightarrow{x})⟩=\sum _j{m}_{j}W(\Vert \overrightarrow{x}-{\overrightarrow{x}}_j\Vert ,h)$$

You can use this cubic spline kernel:

$$W(r,h)=\frac{40}{7\pi h^2}\{\begin{array}{ll}6(q^3-q^2)+1& \text{for }0\le q<\frac{1}{2}\\ 2(1-q{)}^3& \text{for }\frac{1}{2}\le q<1\\ 0& \text{otherwise}\end{array}$$

These equations match what you implemented in problem set 3.

2.2 Pressure Computation

In WCSPH, pressure for particle $i$ is directly computed from density:

$$p_i=B((\frac{{\rho }_i}{{\rho }_0}{)}^{\gamma }-1)$$

, where ${\rho }_0$ is the default density, $\gamma$ and $B$ are constants. When $\gamma =1$, it is actually the pressure for ideal gas. To simulate fluids with high resistance to compresssion (like water), we use $\gamma$ = 7 in our assignment.

2.3 Acceleration Computation

From the Navier-Stokes equations for a fluid without viscosity:

$$\frac{D\overrightarrow{u}}{Dt}+\frac{1}{\rho }\mathrm{\nabla }p=\overrightarrow{g}$$

To calculate acceleration and update velocity, we need the pressure gradient. For particle $i$, this is:

$$⟨\mathrm{\nabla }{p}_i⟩={\rho }_i\sum _j{m}_j(\frac{p_i}{{\rho }_i^2}+\frac{p_j}{{\rho }_j^2}){\mathrm{\nabla }}_i{W}_{i,j}$$

, where ${\mathrm{\nabla }}_i{W}_{i,j}=\frac{\partial W(\Vert {\overrightarrow{x}}_i-{\overrightarrow{x}}_j\Vert ,h)}{\partial {\overrightarrow{x}}_i}$. (You might notice that the gradient of pressure here is different from $\mathrm{\nabla }⟨A⟩$ in problem set 3. See here if you want to know more.)

Substituting back into the N-S equation gives the acceleration:

$$\frac{{\text{d}}^2{x}_i}{\text{d}{t}^2}=-\sum _j{m}_j(\frac{p_i}{{\rho }_i^2}+\frac{p_j}{{\rho }_j^2}){\mathrm{\nabla }}_i{W}_{i,j}+g$$

2.4 Boundary Conditions

Boundary conditions are always one of the most annoying interesting parts of fluid simulation. We implement boundaries using fixed boundary particles with density ${\rho }_0$ and pressure $0$. These particles contribute to pressure and density calculations for fluid particles, pushing them away from boundaries. Three layers of boundary particles ensure that the full kernel support will be covered near boundaries.

2.5 Artificial Viscosity

To improve numerical stability, we add artificial viscosity:

$$a_{viscosity}^i=\{\begin{array}{ll}-\sum _j{m}_j{\mathrm{\Pi }}_{i,j}{\mathrm{\nabla }}_i{W}_{i,j}& v_{i,j}^T{x}_{i,j}<0\\ 0& v_{i,j}^T{x}_{i,j}\ge 0\end{array}$$

Where $v=\frac{\text{d}x}{\text{d}t}$ is velcocity, $v_{i,j}=v_i-v_j$, $x_{i,j}=x_i-x_j$, and ${\mathrm{\Pi }}_{i,j}$ is given as:

$${\mathrm{\Pi }}_{i,j}=-\frac{2\alpha h{c}_s}{{\rho }_i+{\rho }_j}\frac{v_{i,j}^T{x}_{i,j}}{\Vert x_{i,j}{\Vert }^2+ϵh}$$

, where $\alpha$ is the viscosity constant (we use $0.8$ in our code). $ϵ=0.01$ is introduced to avoid singularity when $\Vert x_{i,j}\Vert$ is too small.

The final acceleration for particle $i$ is:

$$\frac{{\text{d}}^2{x}_i}{\text{d}{t}^2}=-\sum _j{m}_j(\frac{p_i}{{\rho }_i^2}+\frac{p_j}{{\rho }_j^2}){\mathrm{\nabla }}_i{W}_{i,j}+g+a_{viscosity}$$

3 Interactive Controls

• W, A, S, D: Control gravity direction
• E, Q: Rotate barrier particles clockwise/counter-clockwise

4 Debugging Tips

• Make sure density and pressure are computed correctly. Density at each particle should be close to ${\rho }_0$, and pressure should not be too large.
• Particles are very energetic without viscosity (as shown here). If you want to see some results or debug before implementing viscosity, consider use a smaller constant $B$ in pressure computation. (a too small $B$ might cause boundary penetration, though)
• If particles stick to boundaries or deform in the first few frames, try setting ${\rho }_i=max({\rho }_0,{\rho }_i)$

5 Submission

You need to include the following in your submission:

1. A PDF file including a link to the chosen commit for your submission. If you submit a link just to your repository, we will assume that you would like us to look at the latest commit on the main branch.
2. A demo video demonstrating the functionality of your simulation.

6 References

• Weakly compressible SPH for free surface flows
• Smoothed Particle Hydrodynamics Techniques for the Physics Based Simulation of Fluids and Solids