FLUID SIMULATION

1
FLUID SIMULATION

• Robert Bridson, UBC
• Matthias Müller-Fischer, AGEIA Inc.

2
Course Details

• www.cs.ubc.ca/rbridson/fluidsimulation
• Schedule
• 145 - 330 Basics Bridson
• 330 - 345 Break
• 345 - 445 Real-Time M-F
• 445 - 530 Extras Bridson

3
The Basic Equations
4
Symbols

• velocity with components (u,v,w)
• ? fluid density
• p pressure
• acceleration due to gravity or animator
• ? dynamic viscosity

5
The Equations

• Incompressible Navier-Stokes
• Momentum Equation
• Incompressibility condition

6
The Momentum Equation
7
The Momentum Equation

• Just a specialized version of Fma
• Lets build it up intuitively
• Imagine modeling a fluid with a bunch of
  particles(e.g. blobs of water)
• A blob has a mass m, a volume V, velocity u
• Well write the acceleration as Du/Dt(material
  derivative)
• What forces act on the blob?

8
Forces on Fluids

• Gravity mg
• or other body forces designed by animator
• And a blob of fluid also exerts contact forces on
  its neighbouring blobs

9
Pressure

• The normal contact force is pressure
  (force/area)
• How blobs push against each other, and how they
  stick together
• If pressure is equal in every direction, net
  force is zeroImportant quantity is pressure
  gradient
• We integrate over the blobs volume(equivalent
  to integrating -pn over blobs surface)
• What is the pressure? Coming soon

10
Viscosity

• Think of it as frictional part of contact
  forcea sticky (viscous) fluid blob resists
  other blobs moving past it
• For now, simple model is that we want velocity to
  be blurred/diffused/
• Blurring in PDE form comes from the Laplacian

11
The Continuum Limit (1)

• Model the world as a continuum
• particles ??Mass and volume ?0
• We want Fma to be more than 00
• Divide by mass

12
The Continuum Limit (2)

• The fluid density is ?m/V
• This is almost the same as the stanford eqn(in
  fact, the form we mostly use!)
• The only weird thing is Du/Dt

13
Lagrangian vs. Eulerian

• Lagrangian viewpoint
• Treat the world like a particle system
• Label each speck of matter, track where it goes
  (how fast, acceleration, etc.)
• Eulerian viewpoint
• Fix your point in space
• Measure stuff as it flows past
• Think of measuring the temperature
• Lagrangian in a balloon, floating with the wind
• Eulerian on the ground, wind blows past

14
The Material Derivative (1)

• We have fluid moving in a velocity field u
• It possesses some quality q
• At an instant in time t and a position in space
  x, the fluid at x has q(x,t)
• q(x,t) is an Eulerian field
• How fast is that blob of fluids q changing?
• A Lagrangian question
• Answer the material derivative Dq/Dt

15
The Material Derivative (2)

• It all boils down to the chain rule
• We usually rearrange it

16
Turning Dq/Dt Around

• For a thought experiment, turn it around
• That is, how fast q is changing at a fixed point
  in space (?q/?t) comes from two things
• How fast q is changing for the blob of fluid at x
• How fast fluid with different values of q is
  flowing past

17
Writing D/Dt Out

• We can explicitly write it out from components

18
D/Dt For Vector Fields

• Say our fluid has a colour variable (RGB vector)
  C
• We still write
• The dot-product and gradient confusing?
• Just do it component-by-component

19
Du/Dt

• This holds even if the vector field is velocity
  itself
• Nothing different about this, just that the fluid
  blobs are moving at the velocity theyre carrying.

20
The Incompressibility Condition
21
Compressibility

• Real fluids are compressible
• Shock waves, acoustic waves, pistons
• Note liquids change their volume as well as
  gases, otherwise there would be no sound
  underwater
• But this is nearly irrelevant for animation
• Shocks move too fast to normally be
  seen(easier/better to hack in their effects)
• Acoustic waves usually have little effect on
  visible fluid motion
• Pistons are boring

22
Incompressibility

• Rather than having to simulate acoustic and shock
  waves, eliminate them from our modelassume
  fluid is incompressible
• Turn stiff system into a constraint, just like
  rigid bodies!
• If you fix your eyes on any volume of space,
  volume of fluid in volume of fluid out

23
Divergence

• Lets use the divergence theorem
• So for any region, the integral of is
  zero
• Thus its zero everywhere

24
Pressure

• Pressure pwhatever it takes to make the
  velocity field divergence free
• If you know constrained dynamics,is a
  constraint, and pressure is the matching Lagrange
  multiplier
• Our simulator will follow this approach
• solve for a pressure that makes our fluid
  incompressible at each time step.

25
Alternative

• To avoid having to solve linear systems, can turn
  this into a soft-constraint
• Track density, make pressure proportional to
  density variation
• Called artificial compressibility
• However, easier to accelerate a linear solver if
  the linear system is well-posed

26
Inviscid Fluids
27
Dropping Viscosity

• In most scenarios, viscosity term is much smaller
• Convenient to drop it from the equations
• Zero viscosity inviscid
• Inviscid Navier-Stokes Euler equations
• Numerical simulation typically makes errors that
  resemble physical viscosity, so we have the
  visual effect of it anyway
• Called numerical dissipation
• For animation often numerical dissipation is
  larger than the true physical viscosity!

28
Aside A Few Figures

• Dynamic viscosity of air
• Density of air
• Dynamic viscosity of water
• Density of water
• The ratio, ?/? (kinematic viscosity) is whats
  important for the motion of the fluid air is
  10 times more viscous than water!

29
The inviscid equations

• A.k.a. the incompressible Euler equations

30
Boundary Conditions
31
Boundary Conditions

• We know whats going on inside the fluid what
  about at the surface?
• Three types of surface
• Solid wall fluid is adjacent to a solid body
• Free surface fluid is adjacent to nothing(e.g.
  water is so much denser than air, might as well
  forget about the air)
• Other fluid possibly discontinuous jump in
  quantities (density, )

32
Solid Wall Boundaries

• No fluid can enter or come out of a solid wall
• For common case of usolid0
• Sometimes called the no-stick condition, since
  we let fluid slip past tangentially
• For viscous fluids, can additionally impose
  no-slip condition

33
Free Surface

• Neglecting the other fluid, we model it simply as
  pressureconstant
• Since only pressure gradient is important, we can
  choose the constant to be zero
• If surface tension is important (not covered
  today), pressure is instead related to mean
  curvature of surface

34
Multiple Fluids

• At fluid-fluid boundaries, the trick is to
  determine jump conditions
• For a fluid quantity q, qq1-q2
• Density jump ? is known
• Normal velocity jump must be zero
• For inviscid flow, tangential velocities may be
  unrelated (jump is unknown)
• With no surface tension, pressure jump p0

35
Numerical Simulation Overview
36
Splitting

• We have lots of terms in the momentum equation a
  pain to handle them all simultaneously
• Instead we split up the equation into its terms,
  and integrate them one after the other
• Makes for easier software design tooa separate
  solution module for each term
• First order accurate in time
• Can be made more accurate, not covered today.

37
A Splitting Example

• Say we have a differential equation
• And we can solve the component parts
• SolveF(q,?t) solves dq/dtf(q) for time ?t
• SolveG(q,?t) solves dq/dtg(q) for time ?t
• Put them together to solve the full thing
• q SolveF(qn, ?t)
• qn1 SolveG(q, ?t)

38
Does it Work?

• Up to O(?t)

39
Splitting Momentum

• We have three terms
• First term advection
• Move the fluid through its velocity field
  (Du/Dt0)
• Second term gravity
• Final term pressure update
• How well make the fluid incompressible

40
Space

• Thats our general strategy in time what about
  space?
• Well begin with a fixed Eulerian grid
• Trivial to set up
• Easy to approximate spatial derivatives
• Particularly good for the effect of pressure
• Disadvantage advection doesnt work so well
• Later particle methods that fix this

41
A Simple Grid

• We could put all our fluid variables at the nodes
  of a regular grid
• But this causes some major problems
• In 1D incompressibility means
• Approximate at a grid point
• Note the velocity at the grid point isnt
  involved!

42
A Simple Grid Disaster

• The only solutions to are
  uconstant
• But our numerical version has other solutions

u
x
43
Staggered Grids

• Problem is solved if we dont skip over grid
  points
• To make it unbiased, we stagger the gridput
  velocities halfway between grid points
• In 1D, we estimate divergence at a grid point as
• Problem solved!

44
The MAC Grid

• From the Marker-and-Cell (MAC) method
  HarlowWelch65
• A particular staggering of variables in 2D/3D
  that works well for incompressible fluids
• Grid cell (i,j,k) has pressure pi,j,k at its
  center
• x-part of velocity ui1/2,jk in middle of x-face
  between grid cells (i,j,k) and (i1,j,k)
• y-part of velocity vi,j1/2,k in middle of y-face
• z-part of velocity wi,j,k1/2 in middle of z-face

45
MAC Grid in 2D
46
Array storage

• Then for a nx?ny?nz grid, we store them as3D
  arrays
• pnx, ny, nz
• unx1, ny, nz
• vnx, ny1, nz
• wnx, ny, nz1
• And translate indices in code, e.g.

47
The downside

• Not having all quantities at the same spot makes
  some algorithms a pain
• Even interpolation of velocities for pushing
  particles is annoying
• One strategy switch back and forth
  (colocated/staggered) by averaging
• My philosophy avoid averaging as much as
  possible!

48
Advection Algorithms
49
Advecting Quantities

• The goal is to solvethe advection equation
  for any grid quantity q
• in particular, the components of velocity
• Rather than directly approximate spatial term,
  well use the Lagrangian notion of advection
  directly
• Were on an Eulerian grid, though, so the result
  will be called semi-Lagrangian
• Introduced to graphics by Stam99

50
Semi-Lagrangian Advection

• Dq/Dt0 says q doesnt change as you follow the
  fluid.
• So
• We want to know q at each grid point at
  tnew(that is, were interested in xnewxijk)
• So we just need to figure outxold (where
  fluid at xijk came from)andq(xold) (what
  value of q was there before)

51
Finding xold

• We need to trace backwards through the velocity
  field.
• Up to O(?t) we can assume velocity field constant
  over the time step
• The simplest estimate is then
• I.e. tracing through the time-reversed flow with
  one step of Forward Euler
• Other ODE methods can (and should) be used

52
Aside ODE methods

• Chief interesting aspect of fluid motion is
  vortices rotation
• Any ODE integrator we use should handle rotation
  reasonably well
• Forward Euler does not
• Instability, obvious spiral artifacts
• Runge-Kutta 2 is simplest half-decent method

53
Details of xold

• We need for a quantity stored at
  (i,j,k)
• For staggered quantities, say u, we need it at
  the staggered location e.g.
• We can get the velocity there by averaging the
  appropriate staggered components e.g.

54
Which u to Use?

• Note that we only want to advect quantities in an
  incompressible velocity field
• Otherwise the quantity gets compressed(often an
  obvious unphysical artifact)
• For example, when we advect u, v, and w
  themselves, we use the old incompressible values
  stored in the grid
• Do not update as you go!

55
Finding q(xold)

• Odds are when we trace back from a grid point to
  xold we wont land on a grid point
• So we dont have an old value of q there
• Solution interpolate from nearby grid points
• Simplest method bi/trilinear interpolation
• Know your grid be careful to get it right for
  staggered quantities!

56
Boundary Conditions

• What if xold isnt in the fluid? (or a nearest
  grid point were interpolating from is not in the
  fluid?)
• Solution extrapolate from boundary of fluid
• Extrapolate before advection, to all grid points
  in the domain that arent fluid
• ALSO if fluid can move to a grid point that
  isnt fluid now, make sure to do semi-Lagrangian
  advection there
• Use the extrapolated velocity

57
Extrapolating u Into Solids

• Note that the inviscid no-stick boundary
  condition just says
• The tangential component of velocity doesnt have
  to be equal!
• So we need to extrapolate into solids,and cant
  use the solids own velocity(result would be bad
  numerical viscosity)

58
Extrapolation

• One of the trickiest bits of code to get right
• For sanity check, be careful about whats known
  and unknown on the grid
• Replace all unknowns with linear combinations of
  knowns
• More on this later (level sets)

59
Body Forces
60
Integrating Body Forces

• Gravity vector or volumetric animator forces
• Simplest scheme at every grid point just add

61
Making Fluids Incompressible
62
The Continuous Version

• We want to find a pressure p so that the updated
  velocityis divergence freewhile
  respecting the boundary conditions

63
The Poisson Problem

• Plug in the pressure update formula into
  incompressibility
• Turns into a Poisson equation for pressure

64
Discrete Pressure Update

• The discrete pressure update on the MAC grid

65
Discrete Divergence

• The discretized incompressibility condition on
  the new velocity (estimated at i,j,k)

66
Discrete Boundaries

• For now, lets voxelize the geometryeach grid
  cell is either fluid, solid, or empty

S
E
E
E
F
F
S
S
F
F
F
F
S
S
F
S
S
F
67
Whats Active?

• Pressure well only solve for p in Fluid cells
• Velocity anything bordering a Fluid cell is good
• Boundary conditions
• p0 in Empty cells (free surface)
• p? whatever is needed in Solid cells so that the
  pressure update makes(Note normal is to grid
  cell!)(Also note different implied pressures
  for each active face of a solid cell - we dont
  store p there)

68
Example Solid Wall

• (i-1,j,k) is solid, (i,j,k) is fluid
• Normal is n(1,0,0)
• Want
• The pressure update formula is
• So the ghost solid pressure is

S
F
69
Discrete Pressure Equations

• Substitute in pressure update formula to discrete
  divergence
• In each fluid cell (i,j,k) get an equation

70
Simplified

• Collecting terms
• Empty neighbors drop zero pressures
• Solid neighbors e.g. ( i-1, j, k) is solid
• Drop pi-1jk - reduce coeff of pijk by 1 (to 5)
• Replace ui-1/2jk in right-hand side with usolid

71
Linear Equations

• End up with a sparse set of linear equations to
  solve for pressure
• Matrix is symmetric positive (semi-)definite
• In 3D on large grids, direct methods
  unsatisfactory
• Instead use Preconditioned Conjugate Gradient,
  with Incomplete Cholesky preconditioner
• See course notes for full details (pseudo-code)
• Residual is how much divergence there is in un1
• Iterate until satisfied its small enough

72
Voxelization is Suboptimal

• Free surface artifacts
• Waves less than a grid cell high arent seen by
  the fluid solver thus they dont behave right
• Left with strangely textured surface
• Solid wall artifacts
• If boundary not grid-aligned, O(1) error it
  doesnt even go away as you refine!
• Slopes are turned into stairs,water will pool on
  artificial steps.
• More on this later

73
Smoke

• As per Fedkiw et al.01

74
Soot

• Smoke is composed of soot particles suspended in
  air
• For simulation, need to track concentration of
  soot on the grid s(x,t)
• Evolution equation
• Extra term is diffusion (if k?0)
• Simplest implementation Gaussian blur
• Boundary conditions
• At solid wall ds/dn0 (extrapolate from fluid)
• At smoke source sssource
• Note lots of particles may be better for
  rendering

75
Heat

• Usually smoke is moving because the air is hot
• Hot air is less dense pressure vs. gravity ends
  up in buoyancy (hot air rising)
• Need to track temperature T(x,t)
• Evolution equation (neglecting conduction)
• Boundary conditions conduction of heat between
  solids and fluids
• Extrapolate T from fluid, add source term

76
Boussinesq and Buoyancy

• Density variation due to temperature or soot
  concentration is very small
• Use the Boussinesq approximationfix
  densityconstant, but replace gravity in momentum
  equation with a buoyancy force
• Constants ? and ? can be tuned

77
Open Boundaries

• If open air extends beyond the grid, what
  boundary condition?
• After making the Boussinesq approximation,
  simplest thing is to say the open air is a free
  surface (p0)
• We let air blow in or out as it wants
• May want to explicitly zero velocity on far
  boundaries to avoid large currents developing

78
Water
79
Where is the Water?

• Water is often modeled as a fluid with afree
  surface
• The only thing were missing so far is how to
  track where the water is
• Voxelized which cells are fluid vs. empty?
• Better where does air/water interface cut
  through grid lines?

80
Marker Particles

• Simplest approach is to use marker particles
• Sample fluid volume with particles xp
• At each time step, mark grid cells containing any
  marker particles as Fluid, rest as Empty or Solid
• Move particles in incompressible grid velocity
  field (interpolating as needed)

81
Rendering Marker Particles

• Need to wrap a surface around the particles
• E.g. blobbies, or Zhu Bridson 05
• Problem rendered surface has bumpy detail that
  the fluid solver doesnt have access to
• The simulation cant animate that detail properly
  if it cant see it.
• Result looks fine for splashy, noisy surfaces,
  but fails on smooth, glassy surfaces.

82
Improvement

• Sample implicit surface function on simulation
  grid
• Even just union of spheres is reasonable if we
  smooth and render from the grid
• Note make sure that a single particle always
  shows up on grid (e.g. radius ?x) or droplets
  can flicker in and out of existence
• If not you must delete stray particles
• But still not perfect for flat fluid surfaces.

83
Level Sets

• Idea drop the particles all together, and just
  sample the implicit surface function on the grid
• In particular, best behaved type of implicit
  surface function is signed distance
• Surface is exactly where ? 0
• See e.g. the book OsherFedkiw02

84
Surface Capturing

• We need to evolve ? on the grid
• The surface (?0) moves at the velocity of the
  fluid (any fluid particle with ?0 should
  continue with ?0)
• We can use same equation for rest of
  volume(doesnt really matter what we do as long
  as sign stays the same at least in continuous
  variables!)

85
Reinitializing

• One problem is that advecting ? this way doesnt
  preserve signed distance property
• Eventually gets badly distorted
• Causes problems particularly for extrapolation,
  also for accuracy in general
• Reinitialize recompute distance to surface
• Ideally shouldnt change surface at all, just
  values in the volume

86
Computing Signed Distance

• Many algorithms exist
• Simplest and most robust (and accurate enough for
  animation) are geometric in nature
• Our example algorithm is based on a Fast Sweeping
  method Tsai02

87
Set Up

• Start by identifying grid cells where ?old
  changes sign (where the surface is)
• Determine points and/or mesh elements to
  approximate surface
• Set nearby grid values to exact distance to
  closest surface element
• Also store index of closest surface element
• Every other grid point is reset to ? Land
  index of closest surface element UNKNOWN
• We then propagate information out to these points

88
Fast Sweeping

• One sweep loop i,j,k through the grid
• Check neighbors for their closest surface element
• If the closest of those is closer than ?ijk
  then update ?ijk and set its index to that point
• Idea keep sweeping in different directions
• i, j, k
• --i, j, k
• i, --j, k
• --i, --j, k
• i, j, --k
• etc

89
Extrapolation

• We often need fluid values extrapolated outside
  of the fluid
• e.g. velocity must be extrapolated out from free
  surface and into solid
• With signed distance this is much simpler

90
Extrapolation (1)

• Interpolate fluid quantities onto stored surface
  points (generated in redistancing)
• Make sure to only use known quantities, not
  values outside the fluid that the physics doesnt
  touch
• Then use index-of-closest-element to copy values
  everywhere needed

91
Extrapolation (2)

• Gradient of distance points parallel to ray to
  closest surface point
• Phrase extrapolation as a PDE
• Then solve the PDE with finite differences
• Careful upwind scheme needed(want to get
  information from points with smaller distance,
  not bigger)