-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathindex.html
253 lines (200 loc) · 14.8 KB
/
index.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
<head>
<style>
body {
padding: 100px;
width: 1000px;
margin: auto;
text-align: left;
font-weight: 300;
font-family: 'Open Sans', sans-serif;
color: #121212;
}
h1, h2, h3, h4 {
font-family: 'Source Sans Pro', sans-serif;
}
</style>
<title>CS 184 Final Project Proposal</title>
<meta http-equiv="content-type" content="text/html; charset=utf-8" />
<link href="https://fonts.googleapis.com/css?family=Open+Sans|Source+Sans+Pro" rel="stylesheet">
<script type="text/javascript" async
src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_CHTML">
// MathJax.Hub.Config({
// jax: ["input/TeX","output/HTML-CSS"],
// displayAlign: "left"
// });
</script>
</head>
<body>
<h1 align="middle">CS 184: Computer Graphics and Imaging, Spring 2017</h1>
<h1 align="middle">Fluid Simulation Project Report</h1>
<h2 align="middle">Team Members: Hubert Jung, Akshay Madhani</h2>
<br>
<div>
<h2 align="left">Abstract</h2>
<p>
In our project we implement fluid simulation by modeling fluids as particles as described in <a href="http://mmacklin.com/pbf_sig_preprint.pdf">Position Based Fluids</a> (Macklin and Müller). This simulation method uses Position Based Dynamics and is based on Smoothed Particle Hydrodynamics (SPH). The big improvement that this method makes over traditional SPH is the ability to efficiently enforce incompressiblity, or constant fluid density. It adds an artificial pressure term to deal with negative pressures and create a surface tension effect. It also includes vorticity confinement to counteract damping, and XSPH viscosity for more coherent motion. The density calculations as well as particle collisions require us to find neighbors of particles, which we are able to compute efficiently using a grid structure.
</p>
<h2 align="left">Overall Loop</h2>
<p>
Our overall simulation loop follows the structure given in Macklin and Müller. This loop can be roughly divided into the following steps:
<ol>
<li> Apply forces and predict new particle positions </li>
<li> Form the particle grid and find neighbors </li>
<li> Repeat for <i>n</i> iterations:
<ul>
<li> Apply incompressiblity constraints and artificial pressure </li>
<li> Calculate particle collisions and environment collisions </li>
</ul>
</li>
<li> Update velocity and position and apply vorticity confinement and viscosity </li>
</ol>
</p>
<h2 align="left">Applying forces</h2>
<p>
For each particle, we calculate the change in velocity due to forces, and predict the change in position due to velocity. We use two simple physics equations:
$$ v_t = v_{t-1} + a{\Delta}t$$
$$ x_t = x_{t-1} + v{\Delta}t$$
</p>
<h2 align="left">Finding neighbors</h2>
<p>
We use a grid-based method to find neighbors. The simulation space is divided into boxes of width <i>w</i> equal to the diameter of a particle. The box a particle is located in is given by the equation:
$$ p_{box} = floor(\frac{p_{particle}}{w}) $$
We form a unique hash from the coordinates of the box and store the particle in a map keyed at that hash.
</p>
<div align="middle">
<img src="images/grid.png" width="600px"/>
<figcaption>Illustration of particles in a grid. Image from <i>Particle Simulation using CUDA</i> (Simon Green)</figcaption>
</div>
<p>
This grid structure greatly reduces the number of particles we need to consider when forming the set of neighbors. For example, for particle collisions we only need to consider the particles located up to one box away. Thus, we would need to iterate over 27 boxes ([-1, 1] in each dimension).
</p>
<h2 align="left">Collisions</h2>
<p>
For each particle, we test for collisions with its neighboring particles and with the environment. If a particle's predicted location is within \(2*radius\) of a neighbor, we simply move it along the vector connecting the two particles so that it is no longer overlapping with the neighbor. We also check to see if the predicted position is on the opposite side of wall from its initial position. If so, we shorten the distance that the particle travels so that it no longer crosses the wall. We additionally scale this distance traveled by a friction parameter, i.e. we multiplied by \((1 - friction)\).
<br> <br>
Note: while we implemented collisions first, they are calculated after all other changes to particle positions.
</p>
<table style="width:100%">
<tr>
<td>
<img src="gifs/basic.gif" width="480px"/>
<figcaption align="middle">Particles falling and colliding</figcaption>
</td>
<td>
<img src="gifs/basic_random.gif" width="480px"/>
<figcaption align="middle">With randomized starting positions</figcaption>
</td>
</tr>
</table>
<h2 align="left">Incompressibility</h2>
<p>
In order to maximize realism, we want to enforce incompressibility. Another way of saying this is that we want the density \(\rho_i\) at each point in the fluid to be equal to a constant rest density \(\rho_o\). This is captured mathematically with the constraint
$$ C_i(\pmb{\textrm{p}}_1,...,\pmb{\textrm{p}}_n) = \frac{\rho_i}{\rho_0} - 1$$
where \(\pmb{\textrm{p}}_1,...,\pmb{\textrm{p}}_n\) represents the set of the particle and its neighbors. We can calculate \(\rho_i\) with the standard SPH density estimator
$$ \rho_i = \sum_j m_jW(\pmb{\textrm{p}}_i - \pmb{\textrm{p}}_j, h)$$
W(<b>r</b>, h) is a smoothing kernel that is only non-zero when \(0 \leq \lvert{\pmb{r}}\rvert \leq h\). Following <a href="http://matthias-mueller-fischer.ch/publications/sca03.pdf"> this paper</a> (Müller et al.), we use the Poly6 kernel for density and the Spiky kernel for the gradient.
$$ W_{poly6}(\pmb{\textrm{r}}, h) = \frac{315}{64\pi h^9}(h^2-{\lvert \pmb{\textrm{r}} \rvert}^2)^3 \quad 0 \leq {\lvert \pmb{\textrm{r}} \rvert} \leq h$$
<figure align="middle">
<img src="images/poly.png" width="200px"/>
<figcaption>Poly6 kernel and its gradient shown by the thick and thin lines, respectively. Image from Müller et al.</figcaption>
</figure>
Our goal in this step is to find a correction vector \(\Delta \pmb{\textrm{p}}\) that can be applied to the position to satisfy the constraint.
$$ C(\pmb{\textrm{p}} + \Delta \pmb{\textrm{p}}) = 0$$
We use Newton's method, stepping along the gradient of the constraint.
$$ \Delta \pmb{\textrm{p}} \approx \nabla C(\pmb{\textrm{p}}) \lambda $$
$$ C(\pmb{\textrm{p}} + \Delta \pmb{\textrm{p}}) \approx C(\pmb{\textrm{p}}) + \nabla C^T \Delta \pmb{\textrm{p}} = 0 $$
Plugging the first equation into the second and solving for lambda yields:
$$ \lambda_i = -\frac{C_i(\pmb{\textrm{p}}_1,...,\pmb{\textrm{p}}_n)}{\sum_k \lvert \nabla_{p_k} C_i \rvert^2 + \varepsilon}$$
where \(\varepsilon\) is an adjustable regularization paramater. The gradient of the constraint is given by
$$\nabla_{p_k} C_i = \frac{1}{\rho_o} \sum_j \nabla_{p_k} W(\pmb{\textrm{p}}_i - \pmb{\textrm{p}}_j, h) $$
You can see from the figure above that the gradient of the Poly6 kernel approaches 0 as r becomes smaller. This would lead to clumping as particles close together would have a very small gradient value. Thus, we instead use the gradient of the Spiky kernel:
$$ W_{spiky}(\pmb{\textrm{r}}, h) = \frac{15}{\pi h^6}(h-{\lvert \pmb{\textrm{r}} \rvert})^3 \quad 0 \leq {\lvert \pmb{\textrm{r}} \rvert} \leq h$$
$$ \nabla W_{spiky}(\pmb{\textrm{r}}, h) = -\frac{45}{\pi h^6}(h-{\lvert \pmb{\textrm{r}} \rvert})^2\frac{\pmb{\textrm{r}}}{{\lvert \pmb{\textrm{r}} \rvert}} \quad 0 < {\lvert \pmb{\textrm{r}} \rvert} \leq h$$
<figure align="middle">
<img src="images/spiky.png" width="200px"/>
<figcaption>Spiky kernel and its gradient shown by the thick and thin lines, respectively. Image from Müller et al.</figcaption>
</figure>
After calculating \(\lambda_i\), we get that the correction vector \(\Delta \pmb{\textrm{p}}\) is
$$ \Delta \pmb{\textrm{p}}_i = \frac{1}{\rho_0} \sum_j(\lambda_i + \lambda_j) \nabla W(\pmb{\textrm{p}}_i - \pmb{\textrm{p}}_j, h). $$
</p>
<figure align="middle">
<img src="gifs/no_vis_no_scorr.gif" width="480px"/>
<figcaption>Fluid with incompressibility</figcaption>
</figure>
<h2 align="left">Tensile Instability</h2>
<p>
One problem with the incompressiblity equations above is that they can cause clumping. When a particle does not have enough neighbors to satisfy the rest density, they will all clump together in an attempt to increase the density. Thus, <i>SPH without a Tensile Instability</i> (Monaghan) adds an artificial pressure term
$$ s_{corr} = -k\bigg(\frac{W(\pmb{\textrm{p}}_i - \pmb{\textrm{p}}_j, h)}{W(\Delta \pmb{\textrm{q}}, h)} \bigg)^n$$
where \(\Delta \pmb{\textrm{q}}\) is some fixed distance within the radius of the smoothing kernel. We played around with parameters and eventually ended with \(k = .001\), \(\lvert \Delta \pmb{\textrm{q}}\rvert = .01h\), and \(n = 4\). Incorporating this term into our \(\Delta \pmb{\textrm{p}}\), our final result is
$$ \Delta \pmb{\textrm{p}}_i = \frac{1}{\rho_0} \sum_j(\lambda_i + \lambda_j + s_{corr}) \nabla W(\pmb{\textrm{p}}_i - \pmb{\textrm{p}}_j, h). $$
In addition to fixing clumping, the artificial pressure also creates a surface tension effect. It causes particle density to be slightly lower than the rest density, which makes particles pull each other closer.
</p>
<figure align="middle">
<img src="gifs/no_vis_with_scorr.gif" width="480px"/>
<figcaption>Fluid with artificial pressure</figcaption>
</figure>
<h2 align="left">Vorticity Confinement and Viscosity</h2>
<p>
Position based dynamics naturally causes damping, and vorticity confinement is a method of adding back in velocity to counteract that damping. Vorticity is related to the rotational energy of the fluid, and can be estimated with the equation
$$ \omega_i = \nabla \times \pmb{\textrm{v}} = \sum_j (\pmb{\textrm{v}}_j - \pmb{\textrm{v}}_i) \times \nabla_{\pmb{\textrm{p}}_j}W(\pmb{\textrm{p}}_i - \pmb{\textrm{p}}_j, h)$$
We can then apply a corrective force
$$ \pmb{\textrm{f}}_i^{vorticity} = \varepsilon\Bigg(\frac{\nabla\lvert\omega_i\rvert}{\lvert\nabla\lvert\omega_i\rvert\rvert} \times \omega_i\Bigg)$$
However, in our project we were unable to implement vorticity confinement correctly. We were unsure of how to calculate the gradient of the magnitude of vorticity, so our final results do not include vorticity confinement at all.
<br>
<br>
On the other hand, we were able to implement viscosity. Essentially, the velocity of each particle is slightly shifted towards the velocities of its neighbors, creating a more coherent and cohesive motion. The equation for viscosity is
$$ \pmb{\textrm{v}}_i^{new} = \pmb{\textrm{v}}_i + c \sum_j (\pmb{\textrm{v}}_j - \pmb{\textrm{v}}_i) * W(\pmb{\textrm{p}}_i - \pmb{\textrm{p}}_j, h)$$
We tested a wide range of values for \(c\), from 0.01 to 0.0001.
<table style="width:100%">
<tr>
<td>
<img src="gifs/high_vis.gif" width="480px"/>
<figcaption align="middle">\(c = .01\)</figcaption>
</td>
<td>
<img src="gifs/low_vis.gif" width="480px"/>
<figcaption align="middle">\(c = .0005\)</figcaption>
</td>
</tr>
</table>
</p>
<h2 align="left">Final Result</h2>
<figure align="middle">
<img src="gifs/4500_particles.gif" width="480px"/>
<figcaption>A full simulation with 4500 particles</figcaption>
</figure>
<h2 align="left">Problems Encountered</h2>
<p>
Besides our failed attempt at implementing vorticity confinement, we faced a variety of problems. This project required a lot of debugging and tinkering with parameters. One common problem we faced early on was particles leaving the box because of NaN positions. Additionally, particles would very easily blow up and bounce all over the box due to a variety of issues. It took very careful adjusting of paramaters to produce good results. Even in our final simulation, there are 2 particles that pop out and bounce around the walls. Another problem we faced is that the particles as a group had a tendency to shift slightly in one direction. Because incompressibility/collision calculations are done in the order of the particle vector, and the vector goes in order from negative in the x direction (and z direction), you can see the entire column of particles shift to the left and to the front. Because we always iterate through the particles in a constant order, all of the particles on one side will shift slightly, and then their neighbors will move to match, and so on. When we first encountered the problem, we decided to move on and return to it; lucky for us, adding the artificial pressure term mostly fixed the problem.
</p>
<h2 align="left">What We Learned</h2>
We learned a lot from this project, including:
<ul>
<li>Kernels: We learned about the purpose of kernels and some of the logic behind their designs. For example, the spiky kernel was created to address the downsides of the poly6 kernel.</li>
<li>Fluids: We learned a lot about fluids and some of the main concerns in making a realistic fluid simulation, such as incompressibility tensile instability.</li>
<li>Fluid simulation is very difficult, even with all the math laid out. We spent a large portion of our time debugging and tweaking paramaters to make the fluid look right. Initially we wanted to finish the simulation early and implement realistic fluid surfaces, but the particle simulation ended up taking all of our time.</li>
</ul>
<h2 align="left">References</h2>
<ul>
<li>Macklin, Miles, and Matthias Müller. "Position based fluids." ACM Transactions on Graphics (TOG) 32.4 (2013): 104.</li>
<li>Müller, Matthias, David Charypar, and Markus Gross. "Particle-based fluid simulation for interactive applications." Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation. Eurographics Association, 2003.</li>
<li>Green, Simon. "Particle Simulation using CUDA." (2012).</li>
<li>Monaghan, Joseph J. "SPH without a tensile instability." Journal of Computational Physics 159.2 (2000): 290-311.</li>
</ul>
<h2 align="left">Video</h2>
<div align="middle">
<iframe width="640" height="360" src="https://www.youtube.com/embed/tLOKEaRWL-o" frameborder="0" allowfullscreen></iframe>
</div>
<h2 align="left">Team Member Contributions</h2>
<p >
Hubert Jung worked primarily on implementing collisions and the grid structure for efficiently finding neighbors. He also helped debug and fix some of the math for incompressibility. Hubert wrote the final website report and created the simulation gifs.
<br> <br>
Akshay Madhani worked primarily on incompressibility, tensile instability, and viscosity. He also experimented with a variety of paramaters to fix problems with particles exploding outwards. Akshay put together the presentation slides and video.
<br> <br>
Special thanks to Chris Correa and Brendan Tang for working with us to determine our constant values.
</p>
</div>
</body>
</html>