Simulating Smoke and Water in CG

1
Simulating Smoke and Water in CG

• Solving the Equations for Fluid Flow.

2
Fundamental Data Representations

• The smoke is represented as a scalar field of
  densities on a grid. The density represents the
  amount of smoke particles in a given voxel
  (typically density is a number between 0 and 1).
• The velocity field is represented as a vector
  field on a grid.
• Both quantities are stored at cell centers.
• The methods we consider are unconditionally
  stable, so we fix the time-step ?t from the
  beginning.

3
Last Time
We talked about how to move smoke (smoke-density)
around
Advection
Diffusion
Sources
Dissipation
4
This Lecture
We will talk about how to make the velocity
evolve in time.
First recall the equations for moving smoke around
The Navier Stokes equations for the velocity
field are given by
Notice the similarity between the two!
5

• Forces

6
Forces I

• Forces accelerate the fluid. In the model we
  consider it is the only way to initiate and
  control the movement of the fluid.

We are solving the equation
Like this
7
Forces II

• Temperature can be represented as a scalar field
  just as the smoke density!
• This means we can have temperature sources,
  advect, diffuse and dissipate temperature and
  just as with smoke density.
• We want temperature to affect the motion of the
  fluid.

8
Forces III

• A combined buoyancy and gravity force looks like
  this

It models that heavy smoke falls to the ground
and hot smoke rises.
? and ? are constants that control the strength
of each term. d is the density of the smoke and T
is the temperature. Tamb is the ambient
temperature. y is the upward direction.
It is up to your imagination to create other
forces. Two examples are the vorticity
confinement force and a wind force consisting of
a main-direction and smaller random
perturbations.
9

• Advection

10
Advecting Fluid Velocity I
Semi-Lagrangian Advection
To find the velocity of a given voxel at
time t?t, we trace the velocity field backwards
in time to time t, and take the velocity from
there.
The path of a particle as a function of time
The Velocity Field
11
Advecting Fluid Velocity II
Semi-Lagrangian Advection
So we will trace the velocity backwards.
We Consider two ways of doing this
p
Let d be any of the components of the fluid
velocity.
1 A first order accurate backwards Euler
time-step
2 A second order accurate backwards Runge Kutta
time-step
To find the velocity of a given voxel at
time t?t, we trace the velocity field backwards
in time to time t, and take the velocity from
there.
12

• Diffusion

13
Density and Velocity Diffusion
The effect of Density Diffusion
Velocity Diffusion makes the fluid move more like
a rigid object.
14
Velocity Diffusion I
Velocity diffusion is modeled by the heat
equation that includes the laplacian. ?? is the
diffusion rate or the viscosity. Let d be one of
the components of the velocity field.
To solve this, we use a finite difference
approximation. To achieve unconditional stability
we need to solve backwards in time. Let the grid
spacing be (?x, ?y, ?z).
15
Velocity Diffusion II
By re-arranging this expression we isolate terms
to time t on the left side and terms to time t?t
on the right side
Where
and
In order to solve this, we must arrange this into
a system of linear equations
Where we order the elements of the grid in
(x,y,z) lexicographic order. We solve for dt?t
and each row in A looks like this
16
Velocity Diffusion III
To solve the linear system of equations
We can observe that the linear system of
equations is sparse, symmetric and positive
definite. This means that the system can be solve
efficiently using Iterative techniques such as
the pre-conditioned conjugate gradient method.
In the project you will be given code for solving
the linear system of equations.
17

• Computing a Divergence Free Velocity Field

18
Computing a Divergence Free Velocity Field I

• We have now solved the equation

But what we want to solve are these equations
Recall the Helmholtz-Hodge Decomposition
Theorem...
19
Helmholtz-Hodge Decomposition Theorem I

• A velocity field V defined on some domain D can
  be uniquely decomposed in the form
• where u is a divergent free velocity field and p
  is the pressure in the fluid.
• Furthermore u is parallel to (the boundary
  of D), i.e. , and the
  boundary condition for the pressure is a Neumann
  boundary condition.

20
Helmholtz-Hodge Decomposition Theorem II
(Slide by Jos Stam, SIGGRAPH 2003)
21
Computing a Divergence Free Velocity Field II
For notational convenience, assume that we set w
equal to u at the beginning of the time-step
Hence we have solved this equation so far
We know from the previous slides that
But how do we find the pressure p?
We apply the divergence on both sides of the
equation...
22
Computing a Divergence Free Velocity Field III
By applying the divergence on both sides of the
equation, we get
Since u is divergence free, this amounts to
If we can solve this equation, we can find the
divergence free velocity field u simply by
subtracting the gradient of the pressure from w
23
Computing a Divergence Free Velocity Field III
So we wish to solve the following equation known
as a Poisson equation
where
And where w is known and p is unknown.
Lets first write out the equation
In this case all the quantities we use are
defined at time t.
24
Computing a Divergence Free Velocity Field IV

• Next we discretize the equation using finite
  differences

25
Computing a Divergence Free Velocity Field V
Now let
Multiplying the equation by ?x2 and re-arranging
terms, we get
Again there is a linear relations ship, so we can
solve a linear system Of equations of the form
Where A is a coefficient matrix, x is the
pressure we are solving for, and b Is the known
divergence of the velocity field.
26
Computing a Divergence Free Velocity Field VI
In forming the linear system of equations we
order all grid points in (x,y,z) lexicographic
order, ie. (1,5,6) comes before (1,6,2).
The (x,y,z)th element of Ax looks like this
This means that the (x,y,z)th row of A looks like
this
The (x,y,z)th element of the vector x is
The (x,y,z)th element of b is
evaluated at (x,y,z)
27
Computing a Divergence Free Velocity Field VII
In practice we solve the following system of
equations which is sparse, symmetric and positive
definite

• Finally, having solved for the pressure, we
• simply subtract the gradient of the pressure
  (calculated using finite differences) from the
  divergent velocity field w to obtain the
  divergence free velocity field u

28
Boundary Conditions and Moving Objects I

• So far we have ignored the fact that some voxels
  may be part of or adjacent to a boundary, or a
  static or moving object immersed in the
  fluid(obstacle).

• The fluid should not enter obstacles.
• The fluid should be pushed by moving objects.

To properly handle this we need to specify
boundary conditions
29
Identifying Fluid and Obstacle Voxels

• In the case of obstacles immersed in the fluid
  In the beginning of an iteration, visit all
  voxels in lexicographic order and determine their
  status as either a FLUID or OBSTACLE voxel.
• In this loop you also assign a linear index to
  each FLUID voxel that will be used in building
  the equation systems.
• Only the OBSTACLE voxels on the boundary between
  fluid and OBSTACLE are used in the equation
  systems, and they are used as boundary conditions
  only.

30
Boundary Conditions and Moving Objects II

• Neumann boundary conditions for the pressure

is the normal of the obstacle.
where
We consider the normal to be axis aligned, since
it is the normal of a voxel.
This affects the way we form the finite
differences in the linear system of equations for
the Poisson equation.
Assume that the normal points in the direction of
the x axis (they are parallel) and that the
boundary is to the left
becomes
31
Boundary Conditions and Moving Objects III

• This means that the row in the matrix A in the
  Poisson equation changes from

To
Which means
32
Boundary Conditions and Moving Objects IV

• We also need boundary conditions for the
  velocity. These are used in the advection step,
  the diffusion step and the divergence calculation
  step.
• In contrast to the case of calculating the
  pressure, we usually explicitly specify the
  values of the velocity at the boundaries. These
  are called Dirichlet Boundary Conditions.

33
Boundary Conditions and Moving Objects V
Assume that the velocity V (u,v). There are two
types of boundary conditions to consider 1) The
velocity normal to the boundary, u, and the
velocity tangential to the boundary, v.
Consider this situation for a static boundary
Static Boundary
Fluid

• The tangential velocity boundary conditions are
  called slip conditions.

34
Boundary Conditions and Moving Objects VI

• For a static obstacle we want the velocity
  component normal to the boundary, u,
• to be zero or point away from the boundary face,
  such that no fluid enters the
• boundary.
• There are three obvious ways to do this
• Set ub u 0.
• Leave u unaffected, and set ub -u. The velocity
  component then interpolates
• to zero at the boundary.
• Set ub 0. If u points into the boundary, set it
  to 0, otherwise leave it unaffected.

35
Boundary Conditions and Moving Objects VII

• For a moving obstacle we want the velocity
  component normal to the boundary, u, to be zero
  at the boundary, or be numerially less than and
  point in the same direction as the velocity of
  the obstacle at the boundary face.
• Again there are three obvious ways to do this
• Set ub u uobstacle, where uobstacle is the
  velocity of the obstacle at the boundary face.
  (an approximation is to just set u ub).
• Leave u unaffected and set ub 2uobstacle u.
  The velocity then interpolates to uobstacle at
  the obstacle boundary face.
• Set ub uobstacle (approximation is to leave ub
  unaffected) and set u uobstacle if
  u-uobstaclelt0 and otherwise leave u unaffected.

Leaving both u and ub unaffected will not
guarantee that fluid will not leak into the
boundary when we are using a Neumann boundary
condition for the pressure.
36
Boundary Conditions and Moving Objects VIII
Tangential boundary conditions for static objects.
Free Slip Fluid flows unaffected
37
Boundary Conditions and Moving Objects IX
Tangential boundary conditions for static objects.
No Slip More turbulent flow along boundaries
38
Boundary Conditions and Moving Objects X
Tangential boundary conditions for static objects.
Halfway between free slip and no slip.
39
Boundary Conditions and Moving Objects XI

• For moving objects the tangential velocity can
  just be set to that of the object or be combined
  with the slip conditions Houston et al.
• The velocity of voxels inside and not on the
  boundary of moving objects should be set to the
  velocity of the moving object.
• You can create effects such as open windows or
  fans by setting the velocity at the boundary
  parallel to the normal of the boundary.
• If you want your fluid to flow out of your
  bounding box just use Neumann boundary conditions.

40
Boundary Conditions and Moving Objects XII

• We can also set boundary conditions for the
  temperature. If the flow should be unaffected,
  use a Neumann boundary condition, otherwise set
  the temperature explicitly to create hot or cold
  objects affecting the smoke.

41
Vorticity Confinement I
(Slide by Jos Stam, SIGGRAPH 2003)
42
Vorticity Confinement II
(Slide by Jos Stam, SIGGRAPH 2003)
43
Vorticity Confinement III
(Slide by Jos Stam, SIGGRAPH 2003)
44
Next Time...

• Vortex Particle Method
• How to write a Research Paper
• Project Description and Handout