Fluid Animation Christopher Batty November 17, 2011 What distinguishes fluids? What distinguishes fluids? • • • • No “preferred” shape Always flows when force is applied Deforms to fit its container Internal forces depend on velocities, not displacements Watch fluid simulation clips Basic Theory Eulerian vs. Lagrangian Lagrangian: Point of reference moves with the material. Eulerian: Point of reference is stationary. e.g. Weather balloon (Lagrangian) vs. weather station on the ground (Eulerian) Eulerian vs. Lagrangian Consider an evolving temperature field. • Lagrangian: Set of moving particles, each with a temperature value. Eulerian vs. Lagrangian Consider an evolving temperature field. • Eulerian: A fixed grid of temperature values, that temperature flows through. Relating Eulerian and Lagrangian Consider the temperature (, ) at a point following a given path, (). (1 ) ( ) (0 ) Two ways temperature changes: 1. Heating/cooling occurs at the current point. 2. Following the path, the point moves to a cooler/warmer location. Time derivatives Mathematically: , = + Chain rule! = + ∙ Definition of = + ∙ Choose = Material Derivative This is called the material derivative, and denoted . Change at a point moving with the velocity field. Change due to movement. = + Change at the current (fixed) point. Advection To track a quantity T simply moving (passively) through a velocity field: =0 or equivalently + = 0 This is the advection equation. Think of colored dye or massless particles drifting around in fluid. Watch dye in fluid clip Equations of Motion • For general materials, we have Newton’s second law: F = ma. • Specialized for fluids: “Navier-Stokes Equations” Navier-Stokes Density × Acceleration = Sum of Forces = Expanding the material derivative… = − ∙ + What are the forces on a fluid? Primarily for now: • Pressure • Viscosity • Simple “external” forces – (e.g. gravity, buoyancy, user forces) Also: • Surface tension • Coriolis • Possibilities for more exotic fluid types: – Elasticity (e.g. silly putty) – Shear thickening / thinning (e.g. “oobleck”, ketchup) – Electromagnetic forces • Sky’s the limit... Watch oobleck clip In full… = − ∙ + Change in velocity at a fixed point Advection Forces (pressure, viscosity gravity,…) Operator splitting Break the equation into steps: 1. 2. 3. 4. Advection: = − ∙ Pressure: = Viscosity: = External: = ℎ 1. Advection Advection Previously, we considered advection of a passive scalar quantity, , by velocity . = − ∙ In Navier-Stokes we saw: = − ∙ Velocity is advected by itself! Advection i.e., The scalar components , of the vector are each advected in the same way. We can reuse the same numerical method. 2. Pressure Pressure What does pressure do? – Enforces incompressibility. Typical fluids (mostly) do not compress. • Exceptions: high velocity, high pressure, … Incompressibility Compressible velocity field Incompressible velocity field Incompressibility Intuitively, net flow into/out of a given region is zero (no sinks/sources). Integrate the flow across the boundary of a closed region: ∙=0 Ω Incompressibility ∙=0 Ω By divergence theorem: ∙=0 But this is true for any region, so ∙ = everywhere. Pressure How does pressure enter in? – Pressure is the force required to ensure that ∙ = 0. – Pressure force has the following form: = − Helmholtz-Hodge Decomposition Input Velocity field Curl-Free (irrotational) Divergence-Free (incompressible) + = = + × = + Aside: Pressure as Lagrange Multiplier Interpret as an optimization: Find the closest to where ∙ = 0 min − s. t. ∙ = 0 A Lagrange multiplier to enforce the constraint is exactly pressure. 3. Viscosity Viscosity What characterizes a viscous liquid? • “Thick”, damped behaviour • Higher resistance to flow Watch viscous honey clip Viscosity Loss of energy due to internal friction between molecules moving at different velocities. Interactions between molecules in different layers causes shear stress that… • acts to oppose relative motion. • causes an exchange of momentum Viscosity Loss of energy due to internal friction between molecules moving at different velocities. Interactions between molecules in different layers causes shear stress that… • acts to oppose relative motion. • causes an exchange of momentum Viscosity Loss of energy due to internal friction between molecules moving at different velocities. Interactions between molecules in different layers causes shear stress that… • acts to oppose relative motion. • causes an exchange of momentum Viscosity Loss of energy due to internal friction between molecules moving at different velocities. Interactions between molecules in different layers causes shear stress that… • acts to oppose relative motion. • causes an exchange of momentum Viscosity Imagine fluid particles with general velocities. Each particle interacts with nearby neighbours, exchanging momentum. Viscosity Amount of momentum exchanged is proportional to: • Velocity gradient, . • Viscosity coefficient, For any closed region, net in/out flow of momentum: ∙ = ∙ Ω So the viscosity force is: = ∙ Viscosity The end result is a smoothing of the velocity. This is exactly the action of the Laplacian operator, ∙ . For scalar quantities, the same operator is used to model diffusion. 4. External Forces External Forces Any other forces you may want. • Simplest is gravity: – = for = (0, −9.81) • Buoyancy models are similar, – e.g., = ( − ) Numerical Methods 1. Advection Advection of a Scalar • First, we’ll consider advecting a scalar quantity, . – temperature, color, smoke density, … • The grid stores and . Eulerian Approximate derivatives with finite differences. + ∙ = 0 FTCS = Forward Time, Centered Space: +1 − +1 − −1 + =0 ∆ 2∆ Unconditionally Unstable! Lax: Conditionally +1 − (+1 + −1 )/2 +1 − −1 + =0 Stable! ∆ 2∆ Many possible methods, stability can be a challenge. Lagrangian • Advect data forward from grid points by integrating position according to grid velocity. ? • Problem: New data position doesn’t necessarily land on a grid point. Semi-Lagrangian • Look backwards in time from grid points, to see where data is coming from. • Interpolate data at previous time. Semi-Lagrangian - Details 1. Determine velocity , at grid point,. 2. Integrate position for a timestep of −∆. e.g. = , − ∆, 3. Interpolate at , call it . 4. Assign , = for the next time. Unconditionally stable! Advection of Velocity • This is great for scalars. What about velocity advection? • Same method: = − ∙ – Trace back with current velocity – Interpolate velocity (component) at that point – Assign it to the grid point at the new time. • Caution: Do not overwrite the velocity field you’re using to trace back! 2. Pressure Recall… Helmholtz-Hodge Decomposition Input Velocity field Curl-Free (irrotational) Divergence-Free (incompressible) + = = + × = + Pressure Projection - Derivation (1) = − Discretize (1) in time… and (2) ∙ = 0 ∆ = − Then plug into (2)… ∆ ∙ − = 0 Pressure Projection Implementation: 1) Solve a linear system for p: ∆ ∙ = ∙ 2) Update grid velocity with: ∆ = − Implementation ∆ ∙ = ∙ Discretize with finite differences: +1 +1 −1 e.g., in 1D: +1 − − −1 − ∆ +1 − ∆ ∆ = ∆ ∆ Solid Boundary Conditions Free Slip: ∙ = 0 52 Staggered vs. Non-staggered • Two choices for grids: – Velocity and pressure co-located +1 −1 +1 −1 – Velocity and pressure staggered +1 +1 −1 Why prefer staggered? • Finite differences line up nicely. – and match, and ∙ match +1 +1 • …which improves stability! −1 3. Viscosity Viscosity PDE: = ∙ Use finite differences. Discretized in time: = + ∆ ∙ ∗ Viscosity – Time Integration = + ∆ ∙ ∗ Explicit integration, let ∗ be : – Just estimate ∆ ∙ – Add to current . – Very unstable. Implicit integration, let ∗ be : – Stable for high viscosities. – Need to solve a system of equations. Viscosity – Implicit Integration Solve for : − ∆ ∙ = Apply separately for each velocity component. e.g. in 1D: +1 − −1 +1 − −1 − ∆ ∆ ∆ − = ∆ Solid Boundary Conditions No-Slip: = 0 Watch no-slip clip 59 4. External Forces Gravity Discretized form is: = + ∆ Just add a downward acceleration to each velocity, scaled by timestep. Gravity • Note: in a closed fluid-filled container, gravity (alone) won’t do anything! – Incompressibility cancels it out. – (Assuming constant density) Start After gravity step After pressure step Simple Buoyancy Look up at grid point, increment velocity based on difference from a reference velocity. = + ∆ − dictates the strength of the buoyancy force. See paper for a more elaborate model: “Visual simulation of smoke”, 2001. User Forces Add whatever additional forces we want, such as: • near the mouse when the user clicks • Paddle forces in Plasma Pong Watch clip Ordering of Steps Order turns out to be important. Why? 1) Incompressibility is not guaranteed at intermediate steps 2) Advecting with a compressible field causes volume/material loss or gain! Ordering of Steps Consider advection in this field: Best ordering 1) Advection 2) Viscosity 3) Add Forces 4) Make incompressible (Pressure) Useful References • Robert Bridson’s SIGGRAPH 2007 Fluid Course notes: http://www.cs.ubc.ca/~rbridson/fluidsimulation/ • “Stable Fluids” [Stam 1999] • “Real-Time Stable Fluid Dynamics for Games” [Stam, 2003] Liquids Liquids What’s missing so far? We need: • A surface representation • Boundary conditions at the surface The Big Picture Velocity Solver Advect Velocities Add Viscosity Add Gravity Project Velocities to be Incompressible What about liquids? Velocity Solver Advect Velocities Add Viscosity Add Gravity Project Velocities to be Incompressible Surface Tracker Surface Tracker Given: liquid surface geometry, velocity field, timestep Compute: new surface geometry by advection. +1 Surface Tracker Ideally: • Efficient • Accurate • Handles merging/splitting (topology changes) • Conserves volume • Retains small features • Smooth surface for rendering • Provides convenient geometric operations • Easy to implement… Very hard (impossible?) to do all of these at once. Surface Tracking Options 1. 2. 3. 4. 5. Particles Level sets Volume-of-fluid (VOF) Triangle meshes Hybrids (many of these) Particles [Zhu & Bridson 2005] Particles Perform passive Lagrangian advection on each particle. For rendering, need to reconstruct a surface. Level sets [Losasso et al. 2004] Level sets Each grid point stores signed distance to the surface (inside <= 0, outside > 0). Surface is interpolated zero isocontour. >0 <= 0 Volume of fluid [Mullen et al 2007] Volume-Of-Fluid Each cell stores fraction f ϵ 0,1 indicating how empty/full it is. Surface is transition region, f ≈ 0.5. 1 1 1 1 1 0. 8 0 0 0 1 0.4 0 0 1 1 0.4 0 0 1 1 0.8 0 0 0 1 1 0.5 0 0 0 Meshes [Brochu et al 2010] Meshes Store a triangle mesh. Advect its vertices, and correct for collisions.