Fluid dynamics Fields Domain R 2 Scalar field f : R Vector - - PowerPoint PPT Presentation

Fluid dynamics

• Domain
• Scalar field
• Vector field

Fields

Ω ⊆ R2 f : Ω → R f : Ω → R2

Types of derivatives

• Derivatives measure how something changes

due to its parameters

• Temporal derivatives
• Spatial derivatives
• gradient operator
• divergence operator
• Laplacian operator

∂f ∂t

∇f = ∂f ∂x, ∂f ∂y T ∇ · f = ∂f x ∂x + ∂f y ∂y ∇2f = ∇ · (∇f) = ∂2f ∂x + ∂2f ∂y

Spatial derivatives

• Gradient: a vector pointing at the steepest uphill
• f the function
• Divergence: measures the net flux
• Laplacian: measures the difference from the

average of neighbors

Quiz

• What is the laplacian of f(x, y) = 2xy +

y2?

• (2y, 2x+2y)
• 2
• x + 2y
• 4

Quiz

What’s the gradient at (π

3 , −π 4 )? What about Laplacian?

Representation

Domain Ω Density ρ : Ω → [0, 1] Velocity u : Ω → R3

Navier-Stoke equations

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f ∇ · u = 0

u: velocity p: pressure s: kinematic viscosity f: body force ρ: fluid density

Momentum equation Mass equation (Incompressibility condition)

Momentum equation

• Each particle represents a little blob of fluid

with mass m, a volume V, and a velocity u

• The acceleration of the particle
• The Newton’s law

a ≡ Du Dt mDu Dt = F

Forces acting on fluids

• Gravity: mg
• Pressure: -∇p
• Imbalance of higher pressure
• Viscosity: µ∇·∇u
• Force that makes particle moves at average

speed

• dynamic viscosity coefficient: µ

Momentum equation

ρDu Dt = ρg − ∇p + µ∇ · ∇u mDu Dt = mg − V ∇p + V µ∇ · ∇u Du Dt + 1 ρ∇p = g + s∇ · ∇u

The movement of a blob of fluid Divide by volume Rearrange equation giving Navier-Stoke

Quiz

• If we want to add an additional force F to the

entire domain of the fluid, how does that change the momentum equation?

Du Dt + 1 ρ∇p = g + s∇ · ∇u

Lagrangian v.s. Eulerian

• Lagrangian point of view describes motion as

points travel around space over time

• Eulerian point of view describes motion as the

change of velocity field in a stationary domain

• Lagrangian approach is conceptually easier, but

Eulerian approach makes the spatial derivatives easier to compute/approximate

Material derivative

Dq Dt = ∂q ∂t + u ∂q ∂x + v ∂q ∂y + w ∂q ∂z = ∂q ∂t + ∇q · u = Dq Dt d dtq(t, x) = ∂q ∂t + ∇q · dx dt

Advection

• Advection describe how quantity, q, moves

with the velocity field u

• Density
• Velocity

Du Dt = ∂u ∂t + u · ∇u Dρ Dt = ∂ρ ∂t + u · ∇ρ

Advection

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Advection Projection Diffusion Body force

Du Dt + 1 ρ∇p = g + s∇ · ∇u

momentum equation

Density advection

∂ρ ∂t = −(u · ∇)ρ

integrate density field using ∂ρ

∂t

compute derivative of density field:

Velocity advection

∂ρ ∂t = −(u · ∇)ρ

∂u ∂t = (u · r)u

integrate velocity field using ∂u

∂t

compute derivative of velocity field: compute derivative of density field: integrate density field

Projection

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Advection

∇ · u = 0

Projection Diffusion Body force

∇ · u > 0? ∇ · u < 0?

Divergence free

∇ ∙ u = 0

Divergence free

∇ ∙ u ≠ 0

Projection

s.t. ∇ · u = 0

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Diffusion

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Advection

∇ · u = 0

Projection Diffusion Body force

High viscosity fluids

Dropping viscosity

• Viscosity plays a minor role in most fluids
• Numeric methods introduce errors which can

be physically reinterpreted as viscosity

• Fluid with no viscosity is called “inviscid”

Body force

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Advection

∇ · u = 0

Projection Diffusion Body force

External forces

• Gravity
• Heat
• Surface tension
• User-specified forces (stirring coffee)

Boundary conditions

• Solid walls
• The fluid cannot go in and out of it
• Control the normal velocity
• Free surfaces
• Air can be represented as a region where

the pressure is zero

• Do not control the velocity at the surface

Solid boundary

• Velocity at the boundary
• No-slip condition:
• Pressure at the boundary
• Make the fluid incompressible AND enforce the

solid boundary condition

u · ˆ n = usolid · ˆ n u = usolid

Physics recap

• Physical quantity represented as fields
• Navier-Stokes PDE describes the dynamics

• What is the representation for fluids?
• use grid structures to represent velocity,

pressure, and other quantities of interest

• What is the equation of motion for fluids?
• approximate Navier-Stokes in a discrete

domain

• compute Navier-Stokes efficiently

Practical challenges

Simple grid structure

pi,j ui,j

(i, j) (i+1, j) (i-1, j) (i, j-1) (i, j+1)

Velocity field

ui,j

ui,j = (ui,j, vi,j)

vector field: u

ui,j

scalar field: u

vi,j

scalar field: v

Computing divergence

ui,j vi,j r · ui,j = ui,j+1 ui,j−1 + vi+1,j vi−1,j 2∆x ui,j+1 ui,j−1 vi−1,j vi+1,j

cell width: ∆x horizontal scalar field vertical scalar field

Derivative evaluation

• Evaluating derivatives using scalar fields

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

∂u ∂t = (u · r)u 1 ρ(rp)x + sr2u + fx ∂v ∂t = (u · r)v 1 ρ(rp)y + sr2v + fy

MAC grid structure

pi,j

pi−1,j

pi+1,j pi,j+1 pi,j−1 vi,j+1/2 vi,j−1/2 ui−1/2,j ui+1/2,j

Quiz

vi,j+1/2 vi,j−1/2 ui−1/2,j ui+1/2,j

How do you compute the divergence at ui,j using MAC grid?

Explicit integration

• Solving for the differential equation explicitly,

namely

• You get this...

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f ut+1 = ut + h ˙ u(t)

Explicit integration

Stable fluids

Invented by Jos Stam Simple, fast, and unconditionally stable

Splitting methods

• Suppose we had a system
• We define simulation function Sf
• Then we could define

∂x ∂t = f(x) = g(x) + h(x) Sf(x, ∆t) : x(t) → x(t) + ∆tf(x) Sf(x, ∆t) : x(t + ∆t) = Sg(x, ∆t) ◦ Sh(x, ∆t)

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Splitting methods

w0(x) w1(x) w2(x) w3(x) w4(x) add force Advect Diffuse Project w0 = u(x, t) u(x, t+Δt) = w4

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Splitting methods

w0(x) w1(x) w2(x) w3(x) w4(x) add force Advect Diffuse Project w0 = u(x, t) u(x, t+Δt) = w4

Body forces

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Splitting methods

w0(x) w1(x) w2(x) w3(x) w4(x) add force Advect Diffuse Project w0 = u(x, t) u(x, t+Δt) = w4

Semi-Lagrangian Advection

Numerical dissipation

• Semi-Lagrangian advection tend to smooth out

sharp features by averaging the velocity field

• The numerical errors result in a different

advection equation solved by semi-Largrangian

• Smooth out small vortices in inviscid fluids

∂q ∂t + u ∂q ∂x = u∆x ∂2q ∂x2

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Splitting methods

w0(x) w1(x) w2(x) w3(x) w4(x) add force Advect Diffuse Project w0 = u(x, t) u(x, t+Δt) = w4

Diffusion

• Solve for the effect of viscosity
• Use an implicit method for stable result

∂w2 ∂t = ν∇2w2 (I − ν∆t∇2)w3(x) = w2(x)

∂u ∂t = −(u · ∇)u − 1 ρ∇p + s∇2u + f

Splitting methods

w0(x) w1(x) w2(x) w3(x) w4(x) add force Advect Diffuse Project w0 = u(x, t) u(x, t+Δt) = w4

Projection

P( )

Projection

• Projection step subtracts off the pressure from

the intermediate velocity field w3

• project(∆t, u) must satisfies two conditions:
• divergence free:
• boundary velocity:

∇ · ut+1 = 0 ut+1 · ˆ n = usolid · ˆ n w4 = w3 − ∆t1 ρ∇p

Divergence-free condition

pi,j

pi−1,j

pi+1,j

pi,j+1

pi,j−1

u1 u2 v1 v2 ut+1

− ut+1

∆x + vt+1

− vt+1

∆x = 0

Replace with

ut+1

u2 − ∆t ρ pi+1,j − pi,j ∆x

resulting in ut+1

= ui+1,j + ui,j 2 − ∆t ρ pi+1,j − pi,j ∆x

Interpolating u2 = ui+1,j + ui,j

2

(from )

Divergence-free condition

pi,j

pi−1,j

pi+1,j

pi,j+1

pi,j−1

Similarly, replace other three terms

u1 u2 v1 v2 ut+1

− ut+1

∆x + vt+1

− vt+1

∆x = 0 ut+1

= ui,j + ui−1,j 2 − ∆t ρ pi,j − pi−1,j ∆x vt+1

= vi,j+1 + vi,j 2 − ∆t ρ pi,j+1 − pi,j ∆x vt+1

= vi,j + vi,j−1 2 − ∆t ρ pi,j − pi,j−1 ∆x

with

Divergence-free condition

Result in a discrete Poisson equation

ut+1

− ut+1

∆x + vt+1

− vt+1

∆x = 0

Assuming time step is 1.0 and density is also 1.0

4pi,j − pi+1,j − pi−1,j − pi,j+1 − pi,j−1 = −0.5∆x (ui+1,j − ui−1,j + vi,j+1 − vi,j−1)

Solve a linear system

p =

A

b

4pi,j − pi+1,j − pi−1,j − pi,j+1 − pi,j−1 = −0.5∆x (ui+1,j − ui−1,j + vi,j+1 − vi,j−1)

4 −1 −1 −1 −1 pi,j

Solve a linear system

• Ap = b
• Construct matrix A:
• diagonal entry: 4
• off-diagonal: -1 if is neighboring cell, 0
• therwise
• A is symmetric positive definite
• Use preconditioned conjugate gradient or Gauss-

Seidel relaxation to solve for p

Gauss-Seidel relaxation

Ax = b A = L + U Lx = b − Ux x(k+1) = L−1(b − Ux(k)) x(k+1)

i

= 1 aii (bi − X

j<i

aijx(k+1)

j

− X

j>i

aijx(k)

j )

Decompose A to Solve U L

Boundary of velocity field

0.2 0.4 0.3 0.2 0.1

• 0.4

0.2

• 0.5

0.2

• 0.5

0.0

• 0.7 -0.8 -0.6 -0.3 -0.1
• 0.7 -0.8 -0.6 -0.3 -0.1

0.2 0.4 0.3 0.2 0.1

• 0.2

0.4 0.5 0.5 0.7

• 0.1
• 0.2
• 0.2

0.0 0.1

u

horizontal velocity

Boundary of velocity field

0.2 0.4 0.3 0.2 0.1

• 0.4

0.2

• 0.5

0.2

• 0.5

0.0

• 0.7 -0.8 -0.6 -0.3 -0.1

0.7 0.8 0.6 0.3 0.1

• 0.2 -0.4 -0.3 -0.2 -0.1

0.2

• 0.4
• 0.5
• 0.5
• 0.7

0.1 0.2 0.2 0.0

• 0.1

v

vertical velocity

Boundary of other fields

0.2 0.4 0.3 0.2 0.1

• 0.4

0.2

• 0.5

0.2

• 0.5

0.0

• 0.7 -0.8 -0.6 -0.3 -0.1
• 0.7 -0.8 -0.6 -0.3 -0.1

0.2 0.4 0.3 0.2 0.1 0.2

• 0.4
• 0.5
• 0.5
• 0.7

0.1 0.2 0.2 0.0

• 0.1

p

• Read: Real-Time Fluid Dynamics for

Games by Jos Stam (Important for Project 5)

• Read: Stable fluids by Jos Stam, Siggraph

99

• Read: SIGGRAPH 2006 course notes on

fluid simulation