Fluid simulation
Marek Trtík
PA199
2
Motivation
Other pictures source: [5]
Picturesource:[3]
Picturesource:[2]
3
 Euler approach:
 Fluid is modelled by a vector field, representing the velocity of the fluid.
 Lagrange approach:
 Fluid is modelled by set of particles.
 Smoothed Particle Hydrodynamics:
 Fluid is modelled by set of particles moved via a velocity vector field.
 Hight-field surface approximation:
 Suitable for simulation of only fluid’s surface, e.g., lake or ocean surface.
Outline
4
Euler approach
5
 Assumptions:
 Incompressible fluid:
 Volume of any subregion of the fluid is constant over time.
 Represented by an incompressible constraint.
 Homogeneous fluid:
 The density of fluid is the same and constant in every region of the fluid and over
time.
 Navier-Stokes equations model a fluid:
 Fluid velocity (motion) represented by a vector field 𝒖 𝒙, 𝑡 .
 Fluid pressure represented by a scalar field 𝑝 𝒙, 𝑡 .
 Partial differential equations define changes in the vector field 𝒖 over time.
Fluid Model
6
 The momentum equations (for each coordinate one):
𝜕𝒖
𝜕𝑡
= − 𝒖 ⋅ ∇ 𝒖 −
1
𝜌
∇𝑝 + 𝜈∇2 𝒖 + 𝒈
 The incompressibility constraint:
∇ ⋅ 𝒖 = 0
 Where (let 𝒙 = 𝑥, 𝑦, 𝑧 ⊤
be a position in space and 𝑡 be simulation time):
 𝒖 𝒙, 𝑡 = 𝑢 𝒙, 𝑡 , 𝑣 𝒙, 𝑡 , 𝑤 𝒙, 𝑡
⊤
is the velocity vector field of the fluid. (computed)
 𝑝(𝒙, 𝑡) is a pressure scalar field of the fluid; used to preserve incompressibility. (computed)
 𝜌 is the density of the fluid, e.g., water 𝜌 = 103 𝑘𝑔
𝑚3.
 𝜈 is the viscosity (resistance to deformation) of the fluid, e.g., honey – high viscosity, water – low
viscosity.
 𝒈(𝒙, 𝑡) is the acceleration vector field of forces acting on the fluid, e.g., gravity 𝒈(𝒙, 𝑡) = 0,0, −10 ⊤ 𝑚
𝑠2.
 ⋅ is the dot product.
Navier-Stokes Equations
advection pressure diffusion external
accel.
7
 Operator of spatial partial derivatives: ∇ =
𝜕
𝜕𝑥
,
𝜕
𝜕𝑦
,
𝜕
𝜕𝑧
⊤
.
 Identifies a direction of a maximum increase of a function at a given time.
 Example: ∇𝑝 =
𝜕𝑝
𝜕𝑥
,
𝜕𝑝
𝜕𝑦
,
𝜕𝑝
𝜕𝑧
⊤
.
 Divergence operator: ∇ ⋅ 𝒖 =
𝜕
𝜕𝑥
,
𝜕
𝜕𝑦
,
𝜕
𝜕𝑧
⊤
⋅ 𝑢, 𝑣, 𝑤 ⊤ =
𝜕𝑢
𝜕𝑥
+
𝜕𝑣
𝜕𝑦
+
𝜕𝑤
𝜕𝑧
.
 Can only be applied to a vector field.
Gradient, Divergence and Laplacian
∇ ⋅ 𝒖 < 0 ∇ ⋅ 𝒖 > 0 ∇ ⋅ 𝒖 = 0
8
 Directional derivative: 𝒖 ⋅ ∇ = 𝑢, 𝑣, 𝑤 ⊤ ⋅
𝜕
𝜕𝑥
,
𝜕
𝜕𝑦
,
𝜕
𝜕𝑧
⊤
= 𝑢
𝜕
𝜕𝑥
+ 𝑣
𝜕
𝜕𝑦
+ 𝑤
𝜕
𝜕𝑧
.
 Therefore, 𝒖 ⋅ ∇ 𝒖 = 𝑢
𝜕
𝜕𝑥
+ 𝑣
𝜕
𝜕𝑦
+ 𝑤
𝜕
𝜕𝑧
𝑢, 𝑣, 𝑤 ⊤
=
𝑢
𝜕𝑢
𝜕𝑥
+ 𝑣
𝜕𝑢
𝜕𝑦
+ 𝑤
𝜕𝑢
𝜕𝑧
𝑢
𝜕𝑣
𝜕𝑥
+ 𝑣
𝜕𝑣
𝜕𝑦
+ 𝑤
𝜕𝑣
𝜕𝑧
𝑢
𝜕𝑤
𝜕𝑥
+ 𝑣
𝜕𝑤
𝜕𝑦
+ 𝑤
𝜕𝑤
𝜕𝑧
.
 Laplacian operator: ∇2 = ∇ ⋅ ∇ =
𝜕2
𝜕𝑥2 +
𝜕2
𝜕𝑦2 +
𝜕2
𝜕𝑧2.
 Example: ∇2
𝒖 = ∇ ⋅ ∇ 𝒖 =
𝜕2
𝜕𝑥2 +
𝜕2
𝜕𝑦2 +
𝜕2
𝜕𝑧2 𝑢, 𝑣, 𝑤 ⊤
=
𝜕2 𝑢
𝜕𝑥2 +
𝜕2 𝑢
𝜕𝑦2 +
𝜕2 𝑢
𝜕𝑧2
𝜕2 𝑣
𝜕𝑥2 +
𝜕2 𝑣
𝜕𝑦2 +
𝜕2 𝑣
𝜕𝑧2
𝜕2 𝑤
𝜕𝑥2 +
𝜕2 𝑤
𝜕𝑦2 +
𝜕2 𝑤
𝜕𝑧2
.
Gradient, Divergence and Laplacian
9
 Beside the fluid we also simulate other quantities, e.g., smoke density,
temperature.
 Represent any such quantity 𝑞 as another scalar/vector field.
 Add a related equation, how 𝑞 changes in time:
𝜕𝑞
𝜕𝑡
= − 𝒖 ⋅ ∇ 𝑞 + 𝜈𝛻2 𝑞 + 𝑆
 Observe the similarity with the momentum equation:
 Advection: − 𝒖 ⋅ ∇ 𝑞
 Diffusion: 𝜈𝛻2
𝑞
 We do not have pressure term.
 𝑆 can be used to simulate constant inflow of 𝑞 into the fluid.
=> Methods for solving the momentum equation can be also applied for
𝑞 equation.
Adding Custom Quantities
10
 The fluid can collide with:
 Static solid objects, like walls.
 Freely moveable solid objects, like piece of wood in water.
 Another fluid, like oil stain on water surface. (not covered in this lecture)
 Our goal is to prevent the fluid to flow into the solid objects.
 Let 𝒏 𝒙, 𝑡 , 𝒖 𝑠 𝒙, 𝑡 be the normal and velocity of the solid surface.
 The boundary constraint for:
 Low viscosity fluid: 𝒖 𝒙, 𝑡 ⋅ 𝒏 𝒙, 𝑡 = 𝒖 𝑠 𝒙, 𝑡 ⋅ 𝒏 𝒙, 𝑡
 High viscosity fluid: 𝒖 𝒙, 𝑡 = 𝒖 𝑠 𝒙, 𝑡
 We can use boundary condition to model fluid source and/or sink.
Boundary Conditions
11
 Discretize the space into a regular grid.
For each cell 𝑖, 𝑗, 𝑘 we store:
 Fluid velocity: 𝒖𝑖,𝑗,𝑘
 Pressure: 𝑝𝑖,𝑗,𝑘
 Any other field: 𝑞𝑖,𝑗,𝑘
 Discretize boundary conditions:
 Mark cells filled by solid objects, e.g., walls.
Discretize fields
Δx
Δ𝑦
2D grid
12
 Use finite differences to approximate partial derivatives.
 Examples:
∇𝑝𝑖,𝑗,𝑘 =
𝑝 𝑖+1,𝑗,𝑘 − 𝑝 𝑖−1,𝑗,𝑘
2Δ𝑥
,
𝑝 𝑖,𝑗+1,𝑘 − 𝑝 𝑖,𝑗−1,𝑘
2Δ𝑦
,
𝑝 𝑖,𝑗,𝑘+1 − 𝑝 𝑖,𝑗,𝑘−1
2Δ𝑧
⊤
∇ ⋅ 𝒖𝑖,𝑗,𝑘 =
𝒖 𝑖+1,𝑗,𝑘 − 𝒖 𝑖−1,𝑗,𝑘
2Δ𝑥
+
𝒖 𝑖,𝑗+1,𝑘 − 𝒖 𝑖,𝑗−1,𝑘
2Δ𝑦
+
𝒖 𝑖,𝑗,𝑘+1 − 𝒖 𝑖,𝑗,𝑘−1
2Δ𝑧
∇2 𝑞𝑖,𝑗,𝑘 =
𝑞 𝑖+1,𝑗,𝑘 −2𝑞 𝑖,𝑗,𝑘+ 𝑞 𝑖−1,𝑗,𝑘
Δ𝑥2 +
𝑞 𝑖,𝑗+1,𝑘 −2𝑞 𝑖,𝑗,𝑘+ 𝑞 𝑖,𝑗−1,𝑘
Δ𝑦2 +
𝑞 𝑖,𝑗,𝑘+1 −2𝑞 𝑖,𝑗,𝑘+ 𝑞 𝑖,𝑗,𝑘−1
Δ𝑧2
Discretize derivatives
13
 Method of splitting:
 Solve a complex equation by a sequence numerical integrations.
𝑑𝑞
𝑑𝑡
= 𝑓 𝑞 + 𝑔 𝑞 →
ො𝑞 = 𝑞 𝑡
+ Δ𝑡𝑓(𝑞 𝑡
)
𝑞 𝑡+Δ𝑡
= ො𝑞 + Δ𝑡𝑔(ො𝑞)
 The result is equivalent to a single integration:
𝑞 𝑡+Δ𝑡
= ො𝑞 + Δ𝑡𝑔 ො𝑞
= 𝑞 𝑡
+ Δ𝑡𝑓(𝑞 𝑡
) + Δ𝑡𝑔 𝑞 𝑡
+ Δ𝑡𝑓(𝑞 𝑡
)
= 𝑞 𝑡 + Δ𝑡𝑓(𝑞 𝑡) + Δ𝑡(𝑔(𝑞 𝑡) + 𝒪(Δ𝑡))
= 𝑞 𝑡 + Δ𝑡(𝑓(𝑞 𝑡) + 𝑔(𝑞 𝑡)) + 𝒪(Δ𝑡2)
= 𝑞 𝑡 + Δ𝑡
𝑑𝑞
𝑑𝑡
+ 𝒪(Δ𝑡2)
Solving Equations
14
 We solve the momentum equation using the splitting method:
𝜕𝒖
𝜕𝑡
= − 𝒖 ⋅ ∇ 𝒖 −
1
𝜌
∇𝑝 + 𝜈∇2 𝒖 + 𝒈
 Start in the current state:
𝒘 𝟎 𝒙 = 𝒖 𝒙, 𝑡
 Apply external accelerations 𝒈:
𝒘1 𝒙 = 𝒘0 𝒙 + Δ𝑡𝒈
(forward Euler)
 Apply fluid advection − 𝒖 ⋅ ∇ 𝒖:
𝒘2 𝒙 = 𝒘1 𝐩(𝒙, −Δ𝑡)
(method of characteristics)
Solving Equations
The new velocity at 𝒙 is the velocity that the
particle had a time Δ𝑡 ago at the location
𝐩(𝒙, −Δ𝑡) (going backward in time along 𝐩).
Picturesource:[3]
15
 We solve the momentum equation using the splitting method:
𝜕𝒖
𝜕𝑡
= − 𝒖 ⋅ ∇ 𝒖 −
1
𝜌
∇𝑝 + 𝜈∇2 𝒖 + 𝒈
 Apply fluid viscosity 𝜈∇2
𝒖:
𝒘3 𝒙 = 𝒘2 𝒙 + Δ𝑡𝜈∇2
𝒘3 𝒙
(backward Euler)
 Lastly, we must compute the pressure −
1
𝜌
∇𝑝 acceleration s.t. we
remove divergence from 𝒘3, i.e., to satisfy the incompressibility:
∇ ⋅ 𝒖 = 0
Solving Equations
16
 Helmholtz-Hodge Decomposition: Any vector field 𝒘 can be uniquely
decomposed to a vector field 𝒖 and a scalar field 𝑝 satisfying:
𝒘 = 𝒖 + ∇𝑝
where 𝒖 is a divergence free, i.e., ∇ ⋅ 𝒖 = 0.
 When we apply divergence operator to both sides of the equation:
∇ ⋅ 𝒘 = ∇2 𝑝
we get a Poisson equation.
 Due to discretization, we get a sparse system of linear equations
=> Use, for example, Jacobi method.
 We use the computed pressure field to get the resulting fluid velocity:
𝒖(𝒙, 𝑡 + Δt) = 𝒘3(𝒙) − ∇𝑝
Solving Equations
17
Euler approach
DEMO!
https://paveldogreat.github.io/WebGL-Fluid-Simulation/
http://haxiomic.github.io/projects/webgl-fluid-and-particles/
18
Lagrange approach
19
 The fluid is represented by 𝑛 particles {𝒫0, … , 𝒫𝑛−1}.
 Each particle 𝒫𝑖 is defined by:
 Mass: 𝑚𝑖
 Position vector: 𝒑𝑖
 Velocity vector: 𝒖𝑖
 Total external force: 𝒇𝑖
 Newton’s equations of motion for moving particles:

𝑑𝒑 𝑖
𝑑𝑡
= 𝒖𝑖 (3 equations in 3D space)

𝑑𝒖 𝑖
𝑑𝑡
=
𝒇 𝑖
𝑚 𝑖
(3 equations in 3D space)
Particles Simulation
20
 The attribute 𝒇𝑖 of a particle 𝒫𝑖 is a sum of all forces acting on the particle.
 We usually want Earth’s gravity to act on particles:
 Force of a homogenous field: 𝑚𝑖 𝒈
 Typically: 𝒈 = 0,0, −10 ⊤
 Interaction between particles 𝒫𝑖 and 𝒫𝑗 via
Lennard-Jones force:
 Let 𝑑𝑖,𝑗 = |𝒑𝑖 − 𝒑 𝑗| and 𝒅𝑖,𝑗 =
𝒑 𝑖−𝒑 𝑗
𝑑 𝑖,𝑗
.
 𝑭𝑖,𝑗 =
𝑘1
𝑑 𝑖,𝑗
𝑚 −
𝑘2
𝑑 𝑖,𝑗
𝑛 𝒅𝑖,𝑗, 𝑭𝑗,𝑖 = −𝑭𝑖,𝑗
 where typically 𝑘1 = 𝑘2, 𝑚 = 4 and 𝑛 = 2.
External Forces
|𝑭𝑖,𝑗| for 𝑘1 = 𝑘2 = 10, 𝑚 = 4, 𝑛 = 2.
Plottedathttps://www.desmos.com/calculator
21
Lagrange approach
DEMO!
22
Smoothed Particle Hydrodynamics
23
 Simulate fluid using a set of 𝑛 particles, i.e., Lagrange approach.
 Compute forces acting on the particles by Euler approach. How?
 Smooth properties of particles into continuous fields.
 Use a smoothing kernel 𝑊(𝑥), e.g., poly6:
𝑊 𝑥 =
315
64𝜋𝑑9 ቊ 𝑑2
− 𝑥2 3
if 0 ≤ 𝑥 ≤ 𝑑
0 otherwise
 Let 𝐴 be a property of particle. Then continuous
field 𝐴(𝒙) is:
𝐴 𝒙 = ෍
𝑗=0
𝑛−1
𝑚𝑗
𝐴𝑗
𝜌𝑗
𝑊(|𝒙𝑗 − 𝒙|) .
 Example: 𝜌 𝒙 = σ 𝑗=0
𝑛−1
𝑚𝑗
𝜌 𝑗
𝜌 𝑗
𝑊(|𝒙𝑗 − 𝒙|) = σ 𝑗=0
𝑛−1
𝑚𝑗 𝑊(|𝒙𝑗 − 𝒙|).
Smoothed Particle Hydrodynamics
𝑊 𝑥 for 𝑑 = 1.
24
 With the fields defined we can use momentum and incompressibility
equations:
𝜕𝒖
𝜕𝑡
= − 𝒖 ⋅ ∇ 𝒖 −
1
𝜌
∇𝑝 + 𝜈∇2 𝒖 + 𝒈, ∇ ⋅ 𝒖 = 0.
 We simulate particles => mass is conserved => ∇ ⋅ 𝒖 = 0 is not needed.
 Particles automatically move with the fluid => − 𝒖 ⋅ ∇ 𝒖 is not needed.
 So, we only solve:
𝜕𝒖
𝜕𝑡
= −
1
𝜌
∇𝑝 + 𝜈∇2 𝒖 + 𝒈.
 Recall second Newton’s equation of motion:
𝑑𝒖 𝑖
𝑑𝑡
=
𝒇 𝑖
𝑚 𝑖
.
 Therefore,
𝒇 𝑖
𝑚 𝑖
= −
1
𝜌(𝒙 𝑖)
∇𝑝(𝒙𝑖) + 𝜈∇2
𝒖(𝒙𝑖) + 𝒈.
Smoothed Particle Hydrodynamics
25
 The pressure field 𝑝 can be obtained from density field 𝜌 by law of ideal gas:
 𝑝 𝒙 = 𝑘(𝜌 𝒙 − 𝜌0), where 𝑘 is a gas constant and 𝜌0 is the environment pressure.
 Derivatives of any field 𝐴(𝒙):
∇𝐴 𝒙 = ෍
𝑗=0
𝑛−1
𝑚𝑗
𝐴𝑗
𝜌𝑗
∇𝑊(|𝒙𝑗 − 𝒙|) , ∇2 𝐴 𝒙 = ෍
𝑗=0
𝑛−1
𝑚𝑗
𝐴𝑗
𝜌𝑗
∇2 𝑊(|𝒙𝑗 − 𝒙|)
where ∇𝑊(|𝒙𝑗 − 𝒙|) = 𝑊′(|𝒙𝑗 − 𝒙|)
𝒙 𝑗−𝒙
𝒙 𝑗−𝒙
, ∇2 𝑊(|𝒙𝑗 − 𝒙|) = 𝑊′′(|𝒙𝑗 − 𝒙|) +
2𝑊′(|𝒙 𝑗−𝒙|)
𝒙 𝑗−𝒙
.
 Forces between two particles generated by fields ∇𝑝, ∇2 𝒖 should be
symmetric => we usually modify their computation:
∇𝑝 𝒙𝑖 = ෍
𝑗=0
𝑛−1
𝑚𝑗
𝑝𝑖 + 𝑝𝑗
2𝜌𝑗
∇𝑊(|𝒙𝑗 − 𝒙𝑖|) , ∇2
𝒖 𝒙𝑖 = ෍
𝑗=0
𝑛−1
𝑚𝑗
𝒖𝑗 − 𝒖𝑖
𝜌𝑗
∇2
𝑊 𝒙𝑗 − 𝒙𝑖 .
Smoothed Particle Hydrodynamics
26
Height-field surface approximation
27
 We model a fluid surface by a function ℎ(𝑥, 𝑦, 𝑡).
 At a point (𝑥, 𝑦) in the XY plane and in time 𝑡 the function defines fluid height z =
ℎ(𝑥, 𝑦, 𝑡).
 Change of ℎ in time is given by:
𝜕2ℎ
𝜕𝑡2 = 𝑣2
∇2
ℎ
where 𝑣 is the speed of waves in the fluid.
 How to solve the equation?
 Introduce an auxiliary function 𝑞 =
𝜕ℎ
𝜕𝑡
.
 Rewrite the equation into this system:
𝜕𝑞
𝜕𝑡
= 𝑣2
∇2
ℎ,
𝜕ℎ
𝜕𝑡
= 𝑞.
 Discretize (next slide).
Fluid Surface Model
28
 We discretize the functions ℎ, 𝑞 by 2D arrays:
ℎ(𝑥0 + 𝑖Δ𝑥, 𝑦0 + 𝑗Δ𝑦, 𝑡0 + 𝑘Δ𝑡) => ℎ𝑖,𝑗
𝑘
𝑞(𝑥0 + 𝑖Δ𝑥, 𝑦0 + 𝑗Δ𝑦, 𝑡0 + 𝑘Δ𝑡) => 𝑞𝑖,𝑗
𝑘
where
 𝑖, 𝑗 are indices to the arrays.
 Δ𝑥, Δ𝑦 are distances between grid cells in X,Y axes.
 𝑘 simulation step number.
 𝛥𝑡 simulation time step.
 NOTE: Usually, 𝑥0 = 𝑦0 = 𝑡0 = 0.
 We solve the discretized system numerically, e.g., using forward Euler method:
𝑞𝑖,𝑗
𝑘+1
= 𝑞𝑖,𝑗
𝑘
+ Δ𝑡𝑣2 ℎ𝑖+1,𝑗
𝑘
−2ℎ𝑖,𝑗
𝑘
+ℎ𝑖−1,𝑗
𝑘
Δ𝑥2 +
ℎ𝑖,𝑗+1
𝑘
−2ℎ𝑖,𝑗
𝑘
+ℎ𝑖,𝑗−1
𝑘
Δ𝑦2 ,
ℎ𝑖,𝑗
𝑘+1
= ℎ𝑖,𝑗
𝑘
+ Δ𝑡𝑞𝑖,𝑗
𝑘+1
.
Discretize Model
𝑥0, 𝑦0
Δ𝑥
Δ𝑦
𝑖Δ𝑥
𝑗Δ𝑦
29
Hight-field surface approximation
DEMO!
30
[1] W.J. Laan, S. Green, M. Sainz; Screen Space Fluid Rendering with
Curvature Flow; I3D 2009.
[2] S. Green; Screen Space Fluid Rendering for Games; GDC 2010.
[3] Jos Stam; Stable Fluids; ACM Transactions on Graphics, 2001.
[4] R.Bridson, M.Müller; Fluid simulation; SIGGRAPH 2007 course notes.
[5] GPU Gems 3; Chapter 30: Real-Time Simulation and Rendering of 3D
Fluids.
[6] GPU Gems; Chapter 38: Fast Fluid Dynamics Simulation on the GPU;
https://developer.download.nvidia.com/books/HTML/gpugems/gpuge
ms_ch38.html
[7] C. Johanson; Real-time water rendering; Master thesis, Lund
University,2004.
References