FLUID SIMULATION

1
(No Transcript)
2
FLUID SIMULATION

• Robert Bridson, UBC
• (Ronald Fedkiw, Stanford)
• Eran Guendelman, Stanford
• Matthias Müller-Fischer, AGEIA Inc.

3
Course Details

• Web http//www.cs.ubc.ca/rbridson/fluidsimulatio
  n
• Schedule
• 830 - 1015 Basics Bridson
• 1015 - 1030 Break
• 1030 - 1115 Film Guendelman
• 1115 - 1200 Games M-F
• 1200 - 1215 Wrap-up Bridson

4
The Basic Equations
5
Symbols

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

6
The Equations

• Incompressible Navier-Stokes
• Momentum Equation
• Incompressibility condition

7
The Momentum Equation
8
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?

9
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

10
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

11
Viscosity

• Think of it as the 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

12
The Continuum Limit (1)

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

13
The Continuum Limit (2)

• The fluid density is ?m/V
• This is almost the same as the standard form of
  the equation(and in fact, the form well mostly
  use!)
• The only peculiar thing is Du/Dtwhat does that
  mean?

14
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

15
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

16
The Material Derivative (2)

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

17
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

18
Writing D/Dt Out

• We can explicitly write it out from components

19
D/Dt For Vector Fields

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

20
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.

21
The Incompressibility Condition
22
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
• Waves moves 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

23
Incompressibility

• Rather than having to simulate acoustic and shock
  waves, lets eliminate them from our
  modelassume fluid is incompressible
• How do we express this?
• If you fix your eyes on any volume of space,
  volume of fluid in volume of fluid out

24
Divergence

• Lets use the divergence theorem
• So for any region of space, the integral ofis
  zero
• Thus the integrand must be zero everywhere

25
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.

26
Inviscid Fluids
27
Dropping Viscosity

• This course focus on water and air
• Except at very small length-scales, viscosity
  term is much smaller than the rest
• 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
i stands for incompressible
29
Aside A Few Figures

• Dynamic viscosity of air
• Density of air
• Dynamic viscosity of water
• Density of water
• But the ratio, ?/? (kinematic viscosity) is
  whats important for the motion of the fluid
• It turns out air is effectively 10 times more
  viscous than water but both are nearly inviscid
  in common situations.

30
Boundary Conditions
31
Boundary Conditions

• We know whats going on inside the fluid what
  about at the surface?
• This morning two 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)

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
Numerical Simulation Overview
35
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, but not covered today.

36
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)

37
Does it Work?

• Up to O(?t)

38
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

39
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

40
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!

41
A Simple Grid Disaster

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

u
x
42
Staggered Grids

• But the problem is solved if we dont skip over
  grid points
• To make this unbiased and accurate, we stagger
  the grid, put velocities at the halfway point
  between grid points
• In 1D, we estimate divergence at a grid point as
• Problem solved!

43
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-component of velocity ui1/2,jk in middle of
  x-face between grid cells (i,j,k) and (i1,j,k)
• y-component of velocity vi,j1/2,k in middle of
  y-face
• z-component of velocity wi,j,k1/2 in middle of
  z-face

44
MAC Grid in 2D
45
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.

46
Advection Algorithms
47
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

48
Semi-Lagrangian Advection

• Dq/Dt says that the 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 the
  fluid that ends up at xijk came from)andq(xold)
  -- what value of q was there

49
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 be used of course

50
Details of xold

• We need to know 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.

51
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!

52
Finding q(xold)

• Odds are when we trace back from a grid point to
  xold we will not land on another 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!

53
Boundary Conditions

• What if xold isnt in the fluid? (or a nearest
  grid point were interpolating from is not in the
  fluid?)
• Solutionextrapolate from value at nearest point
  in fluid
• Can do this extrapolation before advection at all
  grid points in the domain that arent fluid
• ALSO if fluid can move to a grid point that
  didnt used to be fluid, make sure to do
  semi-Lagrangian advection there
• Use the extrapolated velocity

54
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 do need to extrapolate into solids, and
  cant simply use the solids own velocity
• BUT if you use the viscous no-slip condition

55
Body Forces
56
Integrating Body Forces

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

57
Making Fluids Incompressible
58
The Continuous Version

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

59
The Poisson Problem

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

60
Discrete Pressure Update

• The discrete pressure update on the MAC grid

61
Discrete Divergence

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

62
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
63
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)

64
Example Solid Wall
S
F

• Say (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

65
Discrete Pressure Equations

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

66
Simplified

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

67
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

68
Voxelization is Bad

• Free surface artifacts
• Only accurate to O(?x)
• 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
• Slopes are turned into stairs,water will pool on
  artificial steps.

69
Free Surface Fix

• Use the Ghost Fluid Method e.g. Gibou et al.02
• Say xijk is in the fluid and xij1k is outside
• p0 happens at fraction ? along the way
• Make linear interpolant go through 0 there
• Guard against divide by ?0
• Net effect add 1/ ? to coefficient of pijk,
  modified pressure update

70
Solid Wall Fix

• Not quite as simple
• Bad idea
• Using accurate normal and interpolation, project
  out normal component of velocity
• Then do voxel pressure solve
• This looks OK for some dynamic splashing, but
  fails to get hydrostatic (nothing moving) case
  right!

71
Aside Inelastic Dynamics

• Two particles collide
• Update
• Inelastic impulse dissipate as much energy as
  possible
• I.e.

72
Non-Closed System

• What if second particle is scripted?(infinite
  mass, but moving)
• Dissipate as much energy as possible still, but
  now include work done on second particle

73
Variational Pressure Solve

• Pressure update is simply fluid particle
  interaction
• Incompressible means no stored energyfully
  inelastic
• Thus it must minimize kinetic energy
  workwith constraint p0 on free surface
• Equivalent to PDE form

74
Solid Wall Fix

• Discretize the integral (kinetic energy plus
  work) first, then minimize discrete least squares
  problem
• Avoids explicitly handling tricky unaligned solid
  wall boundaries
• Bottom line same old pressure solve, but
  coefficients of the matrix and divergence are
  weighted by amount of fluid in a cell
• See course notes for details
• Cheapest estimate
• Say xijk is in fluid and xi-1jk is in solid
• Fraction ? along line segment is solid/fluid
  interface(1- ?) xi-1jk ?xijk
• Then weight is 1-? instead of 1

75
Advantages

• Easy to implement(just modified matrix entries)
• Consistent 1st order accuracy
• Handles hydrostatic case perfectly
• Handles sub-grid geometry plausibly
• Well-posed linear least squares problem with
  correct null-space
• Automatically symmetric positive (semi-)definite

76
Smoke

• As per Fedkiw et al.01

77
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 in
  fluid)
• At smoke source sssource
• Note lots of particles may be better for
  rendering

78
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
• Heat conduction (maybe negligible) and radiation
• Boundary conditions whatever temperature the
  solids are.

79
Boussinesq and Buoyancy

• Density variation due to temperature or soot
  concentration is very small
• Simpler to make the Boussinesq
  approximationfix densityconstant, but replace
  gravity in momentum equation with a buoyancy
  force
• ? and ? are non-negative constants to tune

80
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

81
Water
82
Where is the Water?

• Water is often modeled as a fluid with a free
  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?

83
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)
• E.g. use Modified Euler (Runge-Kutta 2)

84
Rendering Marker Particles

• Need to wrap a surface around the particles
• E.g. metaballs/blobbies
• Problem rendered surface has 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 miserably on smooth, glassy surfaces.

85
Improvement

• Sample implicit surface function on the
  simulation grid
• Even just union of spheres is reasonable if we
  render from the grid
• Use marching-cubes-like approach for estimates
  needed in pressure solve
• Fraction where fluid/air is
• 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
• But still not perfect for flat fluid surfaces.

86
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

87
Surface Capturing

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

88
Reinitializing

• One problem is that advecting ? this way doesnt
  preserve signed distance property
• Eventually gets badly behaved if fluid is
  swirling a lot
• Can reinitialize recompute signed distance to
  current interface
• Shouldnt change where the interface is,at least
  in continuous variables

89
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

90
Set Up

• Start by identifying grid cells where ?old
  changes sign (where the surface is)
• Estimate points on the surface (Marching-Cubes-lik
  e interpolation)
• Initialize ?new for nearby grid points based on
  which of these points is closest (and sign of
  ?old)
• Also store index of closest surface point at
  those grid points
• Every other grid point is reset to ?/-L ( where
  L is larger than the grid) and index of closest
  point UNKNOWN
• We then iteratively improve this guess for ?

91
Fast Sweeping

• One sweep loop i,j,k through the grid
• Check neighbors from earlier in loop for their
  closest points
• 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

92
Extrapolation

• We often need fluid values extrapolated outside
  of the fluid
• e.g. velocity must be extrapolated out from free
  surface and into solid
• If we have closest surface point information on
  grid, this becomes very simpleinterpolate at
  surface points,then exterior grid points copy
  the value from their closest surface point

93
Numerical Dissipation
94
Numerical Dissipation

• Semi-Lagrangian advection (or in fact, just about
  any usable Eulerian method) has a
  flawnumerical dissipation
• Easiest to see in the semi-Lagrangian context
• When we advect a field q(x,t), the new values are
  smoothly interpolated at various points from the
  old values
• That interpolation smoothes the field
• Can show the error behaves as if we did perfect
  advection followed by a step of viscosity /
  diffusion / blurring

95
The Problem

• In the limit ?x?0, this error goes to zero
• Problem we cant or wont take the limit
• Ideally we want a grid with only just enough
  resolution to represent details we care about
• We may be forced to use something even coarser if
  computer resources too limited
• Then we have features on the order of a few grid
  cells
• Numerical dissipation very quickly smoothes them
  away

96
Dissipation Example (1)

• Start with a function nicely sampled on a grid

97
Dissipation Example (2)

• The function moves to the left(perfect
  advection) and is resampled

98
Dissipation Example (3)

• And now we interpolate new sample values, and
  ruin it all!

99
The Symptoms

• For velocity fields
• It looks like fluids are too sticky (molasses) or
  implausible length scale (scale model)
• Swirly turbulent detail quickly decays, left with
  just boring bulk motion
• For smoke concentration
• Smoke diffuses into thin air too fast,nice thin
  features vanish
• For level sets
• Water evaporates into thin air

100
The Cure?

• The perfect cure hasnt been found yet