[tbd] Jos Stam, Stable Fluids (1999)

semi-Lagrangian simulation method

Fourier transform technique only works given periodic boundaries.

sparse linear solver based simulation - work with any boundary conditions

 

assume: fuild is bounded by rectangle

 

 

http://www.dgp.toronto.edu/people/stam/reality/Talks/FluidsTalk/FluidsTalk_files/v3_document.htm

 

1. Navier-Stokes Equation

 

<Evolution of density>

 

<Evolution of velocity>

 

+ velocity should conserve mass

 

두 식은 비슷한 형태를 지녔지만 두번째 식의 경우, 비선형(non-linear) 항이 있어 풀기 어렵다.

 

1) (assume velocity known)

시간에 따라 밀도(density)가 어떻게 변화하는지 나타낸다.

우변의 세 항이 각각 영향을 미치는 요인이 된다.

첫번째 항은 밀도가 속도장(velocity field)을 따르는 것을 나타낸다. (담배연기를 불면 입김 방향을 따르는 것처럼)

두번째 항은 밀도가 kappa의 속도로 퍼져나갈 수 있다는 것을 나타낸다.

마지막 항인 S는 외부의 소스에 의해 증가해야 함을 나타낸다. 이는 user interface로부터 제공된다.

 

 

2. Algorithm

 

discretize the entire space into identical voxels

전체 공간을 복셀로 나눈다. (2D에서 3D로의 extend는 another index / another for loop으로 쉽게 해결 가능)

 

initial grid of density 빈 상태로 시작함

세 단계를 거쳐 그리드를 업데이트함 (세 개의 항) : add source - diffuse - move

1) add the sources to the grid : multiply this grid by the time step and add them to the density grid

2) diffusing densities : D^n -> D^(n+1)

immediate neighbors (최인접복셀) 만 고려. adjacent faces를 통해 density exchange

single face the exchange = density flux in - density flux out = kappa * timestep (neighbor - center) / h^2

플럭스가 높다 = 큰 time steps, large diffusion rate(kappa), small grid spacings(h)

Dⁿ⁺¹_i,j = Dⁿ _i,j + k Dt (Dⁿ _i-1,j + Dⁿ _i+1,j + Dⁿ _i,j-1 + Dⁿ _i,j+1 - 4Dⁿ _i,j)/h²

-> 이웃 셀보다 더 먼 범위로 확산할 경우 unstable해지기 때문에 성립하지 않는다.

-> implicit techniques: 역으로 확산되는 (-t) 양을 찾으면 항상 성립

Ax = b (A can be huge but is spase -> requires fast linear solver)

 

 

3) moving densities: 속도는 안다고 가정. 각 time step마다 velocity field를 따라 density를 움직인다.

transfer only between neighbors

역시나 instability가 발생할 수 있다. (Δt|u| > h)

stable implicit method를 사용할 경우 linear system은 non-symmetrical matrix를 가진다

KEY: particles + grids = density를 파티클 샘플링 시켜 파티클을 이동시킨 후 각 타임스텝마다 다시 그리드 센터로 이동한다.

trace particle backwards in time -> locate the four closest cells to the point -> interpolate the density from these cells -> 2개의 grid 필요 (previous time step / new interpolates values)

interpolation을 통하기 때문에 density 값들이 항상 bounded되어있다. -> unconditionally stable

 

----

 

Computing Velocities (Vector)

 

1st-term: non-linear "velocity should move along itself"

2nd-term: "diffusion" viscous-like fluids (점도 클수록 vel filed smooth) 'n' equations to solve for n-Dimension

3rd-term: force = sources in the density solver

 

1) Moving Velocity

위와 마찬가지로 두 개의 grid 필요. trace backwards

interpolate a new velocity at that location from the neighboring grid cells

"semi-Lagrangian technique"

 

 

질량보존 Conservation of Mass = the fluid should conserve mass

Flow into cell = Flow out of the cell

Hodge decomposition : field = mass conserving + gradient field (slope function / height field / single scalar field)

그라디언트 필드를 구하기 위해 Poisson equation을 풀어야 함.

sparse symmetrical linear system (conjugate gradient, Jacobi)

 

 

3. Summary

 

UpdateVelocity(U1,U0,F,visc,dt)

AddForce(U1,U0,F,dt)

Diffuse(U0,U1,visc,dt)

Move(U1,U0,U0,dt)

ConserveMass(U1,dt)

 

Very easy to code. Only need:

Particle tracer + grid interpolator

Linear solver (FISHPAK or CG)

So in summary here is all you need to write a fluid solver: a linear solver for the diffusion and the mass conservation step and a good particle tracer and grid interpolator for the self-advection step.

 

 

Velocity field

add force

advect

diffuse due to viscous friction

force to mass

 

Scalar (density) field

add source

transport

diffuse

dissipate

 

 

 

 

 

Simple stable fluid solver

 

periodic boundaries

Fourier transform: 각 포인트에 벡터를 지정하는 대신, 각 wave number (k,l) 에 벡터를 지정한다.

-N/2 < k, l < N/2 (N: resolution of the grid)

viscosity: low pass filter which dampens the higher frequencies. Fourier space의 중심에서 멀어질수록 화살표는 작아짐.

velocity: wave number에 항상 수직. 따라서 모든 벡터가 중심의 동심원에 수직한다.

"incompressibility equation"

 

60 lines of C code : vel (u, v) / forces (u0, v0) / visc / dt

add force : force grid * dt + vel field

self-advection : (move velocity) trace~~interpolation~

diffusion & mass conservation : transform our field to the Fourier domain