流体仿真方法总结

本文禁止转载

B站：Heskey0

Fluid Simulation(Sig 2006 Course Notes)

siggraph 2007 course notes(online) - 太傻 - 博客园 (cnblogs.com)

fluids_notes.pdf (ubc.ca)

注：本文使用2006版本，部分章节的详细内容请自行看2007版本（包括Level-Set, PIC, Octree, Solid-Fluid Coupling等内容）

Part①: The Basics

Chapter 1. The Equations of Fluids

Navier-Stokes:

\[\frac{\part \vec u}{\part t}+\vec u\cdot\nabla \vec u+\frac1\rho\nabla p=\vec g+\nu\nabla\cdot\nabla\vec u\\ \nabla\cdot\vec u=0 \]

1.1 Symbols

\(\vec u\): velocity

\(\rho\): density

\(p\): pressure.

\(\vec g\): 重力加速度 \((0,-9.81,0)m/s^2\).

\(\nu=\mu/\rho\): kinematic viscosity.

body forces: forces applied throughout the whole body of fluid

1.2 The Momentum Equation

\[\frac{\part \vec u}{\part t}+\vec u\cdot\nabla \vec u+\frac1\rho\nabla p=\vec g+\nu\nabla\cdot\nabla\vec u \]

The first equation is called the momentum equation. This really is good old Newton's equation \(\vec F=m\vec a\) in disguise.

\[\nabla\cdot\vec u=0 \]

The second equation is called the incompressibility condition.

1.3 Lagrangian and Eulerian Viewpoints

material derivative: (\(q\) might be density, or velocity, or temperature, or many other things.)

\[\frac{d}{dt}q(t,\vec x)=\frac{\part q}{\part t}+\nabla q\cdot\vec u\\ =\frac{Dq}{Dt} \]

advection(also called convection(对流) or transport): how the quantity, or molecules, or particles, move with the velocity field

advection equation: just one that uses the material derivative, at its simplest setting it to zero

\[\frac{Dq}{Dt}=0 \]

1.4 Incompressibility

divergence-free: A vector-field that satisfies the incompressibility condition

1.5 Dropping Viscosity

Euler equation: The Navier-Stokes equations without viscosity (数值耗散就相当于模拟了流体的粘性，所以可以直接丢弃viscosity term)

\[\frac{D\vec u}{Dt}+\frac1\rho\nabla p=\vec g\\ \nabla\cdot\vec u=0 \]

1.6 Boundary Conditions

free surfaces: Dirichlet boundary condition (流体顶部与空气接触的地方) (知道边界外的压力为0)
- \(p=0\) for void grids
solid walls: Neumann boundary condition (流体与容器接触的地方) (不知道墙的压力是多少)
- 知道压强的梯度，反推墙的压强

Chapter 2. Overview of Numerical Simulation

2.2 Splitting the Fluid Equation

\[\begin{aligned} \frac{Dq}{Dt}&=0\quad\textrm{(advection)}\\ \frac{\part \vec u}{\part t}&=\vec g\quad\textrm{(body force)}\\ \frac{\part \vec u}{\part t}+\frac1\rho\nabla p&=0\quad s.t.\quad\nabla\cdot\vec u=0\quad\textrm{(pressure/incompressibility)} \end{aligned} \]

our algorithm:

advection equation: advects quantity \(q\).
body force equation: use Forward Euler
pressure part: make \(\vec u\) divergence-free

Basic fluid algorithm:

start with an initial divergence-free velocity field \(\vec u\).
For time step \(n=0,1,2,...\):
- Determine a good time step \(\Delta t\) to go from time \(t_n\) to time \(t_{n+1}\).
- Set \(\vec u^A=\textrm{advect}(\vec u^n,\Delta t,\vec u ^n)\)。
- Add \(\vec u^B=\vec u^A+\Delta t\vec g\).
- Set \(\vec u^{n+1}=\textrm{project}(\Delta t,\vec u^B)\).

2.4 Grids

staggered grid(Marker and Cell(MAC) grid) + central difference:

\[\frac{\part q}{\part x}\approx\frac{q_{i+1/2}-q_{i-1/2}}{\Delta x} \]

collocated grid + central difference: (unbiased, 精度高)，两个二阶泰勒展开

\[\frac{\part q}{\part x}\approx\frac{q_{i+1}-q_{i-1}}{2\Delta x} \]

Chapter 3. Advection Algorithm

solving the advection equation(PDE): \(Dq/Dt=0\).

\[q^{n+1}=\textrm{advect}(\vec u,\Delta t,q^n) \]

The brute force approach to solving \(Dq/Dt\) for a time step is to simply write out the PDE, e.g. in one dimension:

\[\frac{\part q}{\part t}+u\frac{\part q}{\part x}=0 \]

if we use Forward Euler(不稳定, unstable) for the time derivative and an accurate central difference(无法捕获解的高频部分) for the spatial derivative we get:

\[\frac{q_i^{n+1}-q_i^n}{\Delta t}+u^n_i\frac{q^n_{i+1}-q^n_{i-1}}{2\Delta x}=0 \]

即：

\[q^{n+1}_i=q^n_i-\Delta tu_i^n\frac{q^n_{i+1}-q^n_{i-1}}{2\Delta x} \]

We will instead take a different, simpler and physically-motivated approach called the “semi-Lagrangian” method

semi-Lagrangian:

假设粒子 \(q^{n+1}\) 正好运动到了grid marker(grid的中心点)
上一帧 \(q^n\) 的位置反算出来，然后可以通过双线性插值算出quantity

example: using Forward Euler and linear interpolation for the semi-Lagrangian operations. For grid point \(x_i\), the particle traced back to \(x_{P}=x_i-\Delta tu_i\). Assuming this lies in the interval \([x_j,x_{j+1}]\), and letting \(\alpha=(x_P-x_j)/\Delta x\) be the fraction of the interval the point lands in, the linear interpolation is \(q^n_P=(1-\alpha)q_j^n+\alpha q_{j+1}^n\). So our update is:

\[q_i^{n+1}=(1-\alpha)q_j^n+\alpha q_{j+1}^n \]

3.1 Boundary Conditions

the normal component of the fluid velocity had better be equal to the normal component of the solid’s velocity, but apart from in viscous flows, the tangential component can be completely different. Thus we usually interpolate the fluid velocity at the boundary, and don’t simply take the solid velocity. However, for the particular case of viscous flows (or at least, a viscous fluid-solid interaction), we can indeed take the shortcut of just using the solid’s velocity instead.

3.2 Time Step Size

the semi-Lagrangian approach is "unconditionally stable": no matter how big \(\Delta t\) is, we can never blow up(不会出现数值爆炸).

If we want to run at real-time rates regardless of the accuracy of the simulation, we can pick \(∆t\) equal to the frame time.

It has been suggested[FF01] that an appropriate strategy is to limit \(∆t\) so that the furthest a particle trajectory is traced is five grid cell widths:

\[\Delta t\le\frac{5\Delta x}{u_{max}} \]

3.2.1 The CFL condition

The CFL condition is a simple and very intuitive necessary condition for convergence. ( Convergence means that if you repeat a simulation with smaller and smaller \(∆t\) and \(∆x\), in the limit going to zero, then your numerical solutions approch the exact solution)

The “domain of dependence” for a point is the set of points that have an effect on the value of the solution at that point.

the CFL condition: convergence is only possible in general if, in the limit as \(∆x → 0\) and \(∆t → 0\), the numerical domain of dependence for each point contains the true domain of dependence.

For semi-Lagrangian methods, the CFL condition is automatically satisfied.

3.3 Dissipation

semi-Lagragian advection中，interpolation(averaging) step will smooth out sharp features，a process called "dissipation"

When we use the simple semi-Lagrangian method to try to solve the advection equation without viscosity, our results look like we were simulating a fluid with viscosity. It’s called “numerical dissipation” (or numerical viscosity, or numerical diffusion—they all mean the same thing).

Thus we need techniques to fix this. More accurate interpolation is a partial answer, proposed for example in [FSJ01], but can be expensive and is still usually inadequate. Later chapters will propose some more effective solutions.

Chapter 4. Making Fluids Incompressible(Projection)

The project routine will subtract off the pressure gradient from the intermediate velocity field \(\vec u\):

\[\vec u^{n+1}=\vec u-\Delta t\frac1\rho\nabla p \]

so that the result satisfies incompressibility inside the fluid:

\[\nabla\cdot\vec u^{n+1}=0 \]

and satisfies the solid wall boundary conditions:

\[\vec u^{n+1}\cdot\hat n=\vec u_{solid}\cdot\hat n \]

4.1 The Discrete Pressure Gradient

two-dimention example:

\[u^{n+1}_{i+1/2,j}=u_{i+1/2,j}-\Delta t\frac1\rho\frac{p_{i+1,j}-p_{i,j}}{\Delta x}\\ v^{n+1}_{i,j+1/2}=u_{i,j+1/2}-\Delta t\frac1\rho\frac{p_{i,j+1}-p_{i,j}}{\Delta x}\\ u^{n+1}_{i-1/2,j}=u_{i-1/2,j}-\Delta t\frac1\rho\frac{p_{i,j}-p_{i-1,j}}{\Delta x}\\ v^{n+1}_{i,j-1/2}=v_{i,j-1/2}-\Delta t\frac1\rho\frac{p_{i,j}-p_{i,j-1}}{\Delta x} \]

4.2 The Discrete Divergence

the divergence in two dimension:

\[\nabla\cdot\vec u=\frac{\part u}{\part x}+\frac{\part v}{\part y} \]

central difference + MAC grid (two dimension):

\[(\nabla\cdot\vec u)_{i,j}\approx\frac{u_{i+1/2,j}^{n+1}-u_{i-1/2,j}^{n+1}}{\Delta x}+\frac{v_{i,j+1/2}^{n+1}-v_{i,j-1/2}^{n+1}}{\Delta x}=0 \]

4.3 The Pressure Equations

将4.1中的式子带入4.2中得到：

\[\frac1{\Delta x}\left[\right. \left(u_{i+1/2,j}-\Delta t\frac1\rho\frac{p_{i+1,j}-p_{i,j}}{\Delta x}\right)- \left(u_{i-1/2,j}-\Delta t\frac1\rho\frac{p_{i,j}-p_{i-1,j}}{\Delta x}\right)\\ +\left(v_{i,j+1/2}-\Delta t\frac1\rho\frac{p_{i,j+1}-p_{i,j}}{\Delta x}\right)- \left(v_{i,j-1/2}-\Delta t\frac1\rho\frac{p_{i,j}-p_{i,j-1}}{\Delta x}\right) \left.\right]=0 \]

再化简，得到：

\[\frac{\Delta t}{\rho}\left(\frac{4p_{i,j}-p_{i+1,j}-p_{i,j+1}-p_{i-1,j}-p_{i,j-1}}{\Delta x^2}\right)=-\left(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta x}+\frac{v_{i,j+1/2}-u_{i,j-1/2}}{\Delta x}\right) \]

this equation is numerical approximations to the "Poisson" problem \(-\Delta t/\rho\nabla\cdot\nabla p=-\nabla\cdot\vec u\).

Boundary conditions:

If grid cell \((i, j + 1)\) is an air cell, then we replace \(p_{i,j+1}\) in this equation with zero
If grid cell \((i + 1, j)\) is a solid cell, then we replace \(p_{i+1,j}\) with the value we compute from the boundary condition there.

对应的公式：

\[p_{i,j+1}=0\\ p_{i+1,j}=p_{i,j}+\frac{\rho\Delta x}{\Delta t}(u_{i+1/2,j}-u_{\textrm{solid}}) \]

4.3.1 Putting them in matrix-vector form

\[Ap=d \]

in which:

\(A\): each row of sparse matrix \(A\) corresponds to one equation, i.e. one fluid cell

\(p\): a vector consisting of all pressure unknowns

\(d\): a vector consisting of the divergences in each fluid grid cell

书中的符号：

\(A_{(i,j,k),(i+1,j,k)}\): the coefficient of \(p_{i+1,j,k}\) in the equation for grid cell \((i,j,k)\). (has to be equal to \(A_{(i+1,j,k),(i,j,k)}\))
In two dimensions, we will store three numbers at every grid cell: the diagonal entry \(A_{(i,j),(i,j)}\) , and the entries for the neighbouring cells in the positive directions, \(A_{(i,j),(i+1,j)}\) and \(A_{(i,j),(i,j+1)}\).

4.3.2 The Conjugate Gradient Algorithm

书上介绍了PCG，这里默认读者都会该算法

4.3.3 Incomplete Cholesky

some direct solver:

LU factorization: \(A=LU\)
LDLT factorization: \(A=LDL^T\) (A is symmetric)
Cholesky factorization: \(A=LL^T\) (A is SPD)

IC(0): Performing Incomplete Cholesky only allowing nonzeros in \(L\) where there are nonzeros in \(A\) is called Level Zero

split \(A\) up into its strict lower triangle \(F\) and diagonal \(D\): \(A=F+D+F^T\).
Then the IC(0) factor \(L\) is of the form: \(L=FE^{-1}+E\), where \(E\) is a diagonal matrix. (all we need to compute and store are the diagonal entries of \(L\), and we can infer the others just from \(A\)!)

two dimensions:

\[E_{(i,j)}=\sqrt{A_{(i,j)}-(A_{(i-1,j)}/E_{(i-1,j)})^2-(A_{(i,j-1)}/E_{(i,j-1)})^2} \]

three dimensions:

\[E_{(i,j,k)}=\sqrt{A_{(i,j,k)}-(A_{(i-1,j,k)}/E_{(i-1,j,k)})^2-(A_{(i,j-1,k)}/E_{(i,j-1,k)})^2-(A_{(i,j,k-1)}/E_{(i,j,k-1)})^2} \]

In these equations, we replace terms referring to a non-fluid cell (or cell that lies off the grid) with zero.

4.3.4 Modified Incomplete Cholesky

可以作为Preconditioner

4.4 Projection

\(\vec u^{n+1}=\textrm{project}(\Delta t,\vec u)\):

Calculate the divergence \(d\) (the right-hand side) with modifications at solid wall boundaries
Set the entries of \(A\).
Construct the MIC(0) preconditioner(也可以用MultiGrid)
Solve \(Ap=d\) with MICCG(0), i.e. the PCG algorithm with MIC(0) as preconditioner
Compute the new velocities \(\vec u^{n+1}\) according to the pressure gradient update to \(\vec u\).

why this routine is called project?

projection operator: \(P^2=P\). It turns out that our transformation from \(\vec u\) to \(\vec u^{n+1}\) is indeed a linear projection

4.5 Accurate Curved Boundaries

Part②: Different Types of Fluid

Chapter 5. Smoke

To model the most important effects of smoke, we need two extra fluid variables:

\(T\): temperature
\(s\): concentration

buoyancy force:

\[f_{\textrm{buoy}}=(0,-\alpha s+\beta(T-T_{amb}),0) \]

where:

\(T_{\textrm{amb}}\): ambient background temperature (say \(300K\), or \(270K\) in Canada)

\(\alpha,\beta\): nonnegative constants which can be set by the user to achieve different behaviour.

We typically store \(T\) and \(s\) in the grid cell centers (the same as pressure). To add the buoyancy force to the \(v\) component of velocity in its offset MAC grid location, we thus need to average the \(T\) and \(s\) values at the neighbouring centers, e.g. \(T_{i,j+1/2,k} = (T_{i,j,k} + T_{i,j+1,k})/2\).

just as we dropped the viscosity term in the momentum equation because the unavoidable numerical dissipation mimics it anyhow, we need not explicitly model heat and smoke diffusion: they’ll happen in our numerical methods whether we want them or not. Thus our new evolution equations are simply:

\[\frac{DT}{Dt}=0\quad\frac{Ds}{Dt}=0 \]

When we advect velocity in our simulation, we will also advect temperature and smoke concentration.

Boundary conditions:

temperature field:
- We can extend the temperature field inside objects just by taking the actual temperature of the object (and in this way, we naturally get convection off hot objects)
- outside the domain as the ambient background temperature \(T_{\textrm{amb}}\).
concentration:
- we should extrapolate the smoke concentration into the object, i.e. when asked for smoke inside an object we interpolate from the closest values outside.

The trouble with this model is that, when simulated on a typical coarse grid, the numerical dissipation from the simple semi-Lagrangian advection scheme we use is excessive. The interesting features—small vortices, sharp divides between smoky and clear air—simply vanish far faster than desired (使用粗网格时，advection scheme的数值耗散会很大，小涡会快速消失).

5.1 Vorticity Confinement

To counter the problem of vortices artificially disappearing too fast, we add a “vorticity confinement” force which boosts the strength of vortices in the flow, keeping it lively and interesting for much longer. (额外施加一个vorticity confinement力，推动涡以维持涡的存在)

vorticity(涡量): 流体速度场的旋度，方向沿着流场旋转所围绕的轴.

\[\vec\omega=\nabla\times\vec u \]

涡量为正，则逆时针旋转

first step: find where the vortices are.

构建一个unit vector \(\vec N\) 指向涡的中心：

\[\vec N=\frac{\nabla|\vec\omega|}{||\nabla|\vec\omega|||} \]

second step: 施加力.

涡量与 \(\vec N\) 叉乘，即沿着流场旋转的方向施加力促进旋转：

\[f_{\textrm{conf}}=\epsilon\Delta x(\vec N\times\vec\omega) \]

\(\epsilon\) is a parameter that can be adjusted to control the effect of vorticity confinement.

numerical implementation: (central difference)

计算voticity:

\[\vec\omega_{i,j,k}=( \frac{w_{i,j+1,k}-w_{i,j-1,k}}{2\Delta x}-\frac{v_{i,j,k+1}-v_{i,j,k-1}}{2\Delta x},\\ \frac{u_{i,j,k+1}-u_{i,j,k-1}}{2\Delta x}-\frac{w_{i+1,j,k}-w_{i-1,j,k}}{2\Delta x},\\ \frac{v_{i+1,j,k}-v_{i-1,j,k}}{2\Delta x}-\frac{u_{i,j+1,k}-u_{i,j-1,k}}{2\Delta x} ) \]

计算gradient of \(|\vec\omega|\):

\[\nabla|\vec\omega|_{i,j,k}=(\frac{|\vec\omega|_{i+1,j,k}-|\vec\omega|_{i-1,j,k}}{2\Delta x}\\ \frac{|\vec\omega|_{i,j+1,k}-|\vec\omega|_{i,j-1,k}}{2\Delta x}\\ \frac{|\vec\omega|_{i,j,k+1}-|\vec\omega|_{i,j,k-1}}{2\Delta x} ) \]

计算 \(\vec N\): (防止除0)

\[\vec N_{i,j,k}=\frac{\nabla|\vec\omega|_{i,j,k}}{||\nabla|\vec\omega|_{i,j,k}||+10^{-20}} \]

Vorticity is an extremely important fluid quantity, and by taking the curl of the momentum equation one can derive the “vorticity equation” which shows how vorticity evolves in time. (对原始NS求旋度，得到vorticity equation)

5.2 Sharper Interpolation

前面提到：Using the standard bilinear/trilinear interpolation in semi-Lagrangian advection quickly smears out sharp features in these fields, which for large scale smoke (where true physical diffusion is negligible) looks distinctly wrong. (使用线性插值会导致sharp features丢失)

Catmull-Rom spline: Given values \(q_{i-1},q_i,q_{i+1},q_{i+2}\), we want to estimate \(q\) at point \(x\) between \(x_i\) and \(x_{i+1}\):

\[q(x)\approx q_i+d_ix+(3\Delta q-2d_i-d_{i+1})x^2+(-2\Delta q+d_i+d_{i+1})x^3 \]

in which:

\[d_i=\frac{q_{i+1}-q_{i-1}}{2}\\ d_{i+1}=\frac{q_{i+2}-q_i}{2}\\ \Delta q=q_{i+1}-q_i \]
we set \(d_i\) or \(d_{i+1}\) to zero if it has sign different from \(\Delta q\). Then the interpolant is guaranteed to be monotonic, and we cannot go unstable

Chapter 6. Water and Level Sets

6.1 Marker Particles

marker particles: any grid cell containing a marker particle is marked fluid, and the empty cells are air.

After the particles have been moved, we label each grid cell that contains a marker particle as fluid and the rest as empty (or solid). This information is then used in the other steps of the fluid simulation, such as the pressure projection.

Marker particles are extremely simple and robust, but they do have one glaring flaw: it’s difficult to get accurate information about the water surface out of them. They can be rendered with a blobby implicit surface. Blobbies can give a possible answer again, but it is noisy. This motivates turning to an alternative method for capturing the location of the water: level sets.

6.2 Level Sets

an implicit surface is one defined as the set of zeros of a function. Letting the function be \(\phi(\vec x)\), the surface is:

\[\{\vec x|\phi(\vec x)=0\} \]

对于 \(\phi(\vec x)\):

the interior (water in our case) is the set where \(\phi(\vec x)<0\)

the exterior (air) is the set where \(\phi(\vec x)>0\).

level set: an implicit surface whose function \(\phi\) is just represented with values sampled on a grid. (to be precise, we will assume they are sampled at the grid cell centers)

level sets are perfectly suited for using PDE’s to evolve surfaces.

优点：（相比于marker particles）

Level sets are attractive for water simulation since they can very easily model smooth surfaces.
easily provide information such as whether a point \(\vec x\) is inside the fluid

6.2.1 Signed Distance

computing signed distance to facilitate extrapolating the velocity field.

let \(\hat n\) be the unit length direction to the closest point on the surface

可以知道 \(\phi\) 有以下性质：

\[\nabla\phi\cdot\hat n=1\\ \hat n=\nabla\phi\\ ||\nabla\phi||=1 \]

即：由于 \(\phi\) 是距离函数，x移动多少距离，\(\phi\) 就改变多少

6.2.2 Calculating Signed Distance

Algorithm:

Find the closest points on the surface for the nearby grid points, setting their signed distance values accordingly. Set the other grid points as unknown.
Loop over the unknown grid points \((i, j, k)\) in a chosen order:
- Find all neighbouring grid points who have a known signed distance and closest surface point.
- Find the distances from \(\vec x_{i,j,k}\) to those surface points. If \(\vec x_{i,j,k}\) is closer than a neighbour, mark the neighbour as unknown again.
- Take the minimum of the distances, determine if \((i, j, k)\) is inside or outside, and set its signed distance value accordingly

(1) Finding Closest Points

An alternative for marker particles is to simply initialize cells containing particles to \(−1\) and empty cells to \(+1\); this can give stair-step artifacts, so it’s preferably to smooth this \(±1\) field out with weighted averages before calculating signed distance.

方法：

The simplest approximation in this case (grid values already given) is to estimate the surface points along the lines between a point and its nearest neighbours by finding where the linear interpolant is zero
A more sophisticated attack is to estimate the direction to the closest point with the gradient, evaluated with a finite difference, and then do a line search for a zero of the interpolated implicit function along this direction
A third alternative is to first approximate the surface in each grid cell using marching cubes or variants, and then calculate the exact closest point on the polygonal surface

6.2.3 Boundary Conditions

we compute signed distance in the water and air only relative to the water-air free surface, ignoring solid boundaries. We then extrapolate the resulting \(\phi\) into the solid. If done carefully enough this can enable impressive simulations of surface tension(表面张力) interactions.

A simpler approach is simply to extrapolate \(\phi\) out as a constant, as is done with velocity below (section 6.3), then recompute the \(\phi\) signed distance within the solid. (this requires having a signed distance function for the solid, \(ψ\), which will guide the extrapolation of \(\phi\)) (\(\psi\) 内负外正)

6.2.4 Moving the Level Set

The most critical part of the simulation is actually moving the surface represented by the level set.

问题：

Unfortunately, advecting a signed distance field doesn’t in general preserve the signed distance property. Thus we need to periodically (every time step, or possibly every few time steps) recalculate signed distance as we have discussed above (6.2.2).
the worst problem is numerical dissipation. (ripples, droplets等高频features会小时)
The fundamental issue is that as the level set surface moves, we keep resampling the signed distance function on the fixed grid. (但采样不到高频信号)

To fix this, we thus need to augment our grid representation of the level set. This is the subject of chapter 10, so I will say no more here.

6.3 Velocity Extrapolation

Our fluid simulation loop looks like this:

Make sure the fluid/air signed distance function \(\phi^n\) is signed distance if necessary.
Extrapolate the divergence-free velocity field from the last time step into the air region, to get an extended velocity field \(\vec u^n\).
Add body forces such as gravity and vorticity confinement to get an interim velocity field..
Advect the interim velocity field and \(\phi\) in the divergence-free \(\vec u^n\) velocity field.
Solve for and apply pressure to get a divergence-free velocity \(\vec u^{n+1}\) in the new fluid region.

We now will address the second step, the extrapolation of velocity out into the air.

the directional derivative of the extrapolated velocity should be zero in the direction \(\nabla\phi\):

\[\nabla\phi\cdot\nabla u=0\quad\nabla\phi\cdot\nabla v=0\quad\nabla\phi\cdot\nabla w=0 \]

We’ve thus changed the extrapolation condition into a simple first order linear PDE, then we can apply finite difference where we only use closer neighbours.

6.4 Surface Tension

Surface Tension: Water molecules are more attracted to other water molecules than air molecules, and vice versa

The simplest linear model of surface tension can in fact be phrased in terms of a potential energy equal to a surface tension coefficient \(γ\) times the surface area between the two fluids. (\(γ\) for water and air at normal conditions is approximately \(0.073J/m^2\)) The force seeks to minimize the surface area.

书中有pressure的详细推导过程，这里省略

the pressure at the surface of the water is:

\[p=\gamma\nabla\cdot\nabla\phi \]

It turns out that that \(κ = ∇ · ∇\phi\) is termed the “mean curvature” of the surface, a well-studied geometric quantity that measures how curved a surface is.

水的surface tension coefficient：\(\gamma=0.073J/m^2\).

Chapter 7. Fire

Combustion comes in two flavours: “detonation” and “deflagration”

Detonation occurs in high explosives, where the fuel burns so fast that the fluids expand in a supersonic blast wave. Although explosions are a common staple of visual effects, it generally isn’t a good idea to use accurate physicsbased simulation for the detonation part. First of all, real physical detonations aren’t usually visible (they’re moving at supersonic speeds)—it’s their effects on the surroundings that are visible. It’s far easier to directly animate those effects, and counterintuitively is often more convincing since normal audiences don’t have much knowledge of real blast wave mechanics. The most visible part of Hollywood explosions comes from the billowing smoke produced in combustion, which we already have covered.
Deflagration is combustion at regular speeds without shock waves, giving rise to flames. It’s what we normally think of as fire. Nguyen et al. demonstrated a very effective way to model fire physically in [NFJ02].

Chapter 8. Viscous Fluids

8.1 Stress

net force: 合力\(f=ma\), 什么是Net Force（合力）?

traction: \(\vec t(\vec x,\hat n)\), a vector field of where in space we are measuring it and what the normal to the contact surface there is.

It has units of force per area(单位是力除以表面积): to get the net contact force on a volume of fluid \(Ω\), we integrate the traction over its surface:

\[\vec F=\iint_{\part\Omega}\vec t(\vec x,\hat n) \]

stress tensor: rank-two tensor \(\sigma(\vec x)\), only depends on position, and must be symmetric. (Since the unit normal has no units, the stress tensor is also measured as force per area, just like traction.)

\[\sigma=-p\delta \]

where we use \(\delta\) to mean the identity tensor

\[\vec t=\sigma\hat n \]

8.2 Applying Stress

a general momentum equation for a continuum(elastic solids as well as fluids):

\[\frac{D\vec u}{Dt}=\frac1\rho\vec g+\frac1\rho\nabla\cdot\sigma \]

For general fluid flow, pressure is still a very important quantity, so we will explicitly separate it out from the rest of the stress tensor:

\[\sigma=-p\delta+\tau \]

where \(τ\) is also a symmetric tensor. (for more details in the next subsection)

8.3 Strain Rate and Newtonian Fluids

viscosity: as a region of fluid slips past another, momentum is transferred between them to reduce the difference in velocity

How much the dot-product between two vectors changes is thus a measure of how fast the fluid is deforming. The rate of change of dot-products of vectors in the flow is determined by the symmetric matrix strain rate tensor:

\[D = \frac12(\nabla\vec u+\nabla\vec u^T) \]

We are looking for a symmetric tensor \(τ\) to model stress due to viscosity:

\[\tau=2\mu D \]

where \(\mu\) is the coefficient of dynamic viscosity. Fluid that satisfy this relationship are called Newtonian Fluid. (water and air)

divergence theorem:

\[\iint_{\part\Omega}v\hat n=\iiint_\Omega\nabla\cdot v \]

in which:

\(v\): any vector function

\(\Omega\): domain

\(\part\Omega\): boundary of the domain

\(\hat n\): normal vector on the boundary

本质上：domain内部的散度会相互抵消，只在边界上散度才有值，所以domain三重积分等价于boundary二重积分

Plugging \(\sigma,\tau\) into the momentum equation, we get:

\[\frac{D\vec u}{Dt}+\frac1\rho\nabla p=\frac\mu\rho[\nabla\cdot\nabla\vec u+\nabla(\nabla\cdot\vec u)] \]

If the flow is incompressible, i.e. \(\nabla\cdot\vec u=0\), we end up back at NS equation

8.4 Boundary Conditions

the boundary condition for the stress at the free surface is:

\[\sigma\hat n=-p\hat n+\tau\hat n=0 \]

Note that if the viscous stress \(τ\) is zero, this reduces to \(p = 0\) as before

no-slip boundary condition:

\[\vec u=\vec u_{\textrm{solid}} \]

Chapter 9. Non-Newtonian Fluids: Goop and Sand

9.1 Viscoelastic Material: Goop

Regular fluids retain no memory of their original state, whereas visco-elastic fluids can (partially) “bounce back” to the original state when applied force is removing, behaving more like an elastic solid in this respect. (黏糊糊的流体就想弹性固体一样，会恢复原状)

(TOG 2004) A method for animating viscoelastic fluids.

9.2 Granular Materials: Sand

(TOG 2005) Animating sand as a fluid

Part③: More Algorithms

Chapter 10. The Particle-Level-Set Method

(2002) A hybrid particle level set method for improved interface capturing

(TOG 2002) Animation and rendering of complex water surfaces

Chapter 11. Particle-In-Cell Methods

Games 201有讲

Chapter 12. Octrees

Sparse tiled grid structures that do not store anything in empty regions can help for free surface flow, but to efficiently deal with big differences in resolution, an adaptive approach is required. The simplest such method is using octrees (in two dimensions, quadtrees) instead of regular grids.

(TOG 2004) Simulating water and smoke with an octree data structure

Chapter 13. Solid-Fluid Coupling

siggraph course的2006版中有详细介绍，请自行阅读

Part④: Real-Time Fluids

Chapter 14. Real-Time Simulations

14.1 Real-Time Water

为了能够实时解算水，需要一些简化算法：

Procedural Water: sine waves, FFT ocean
Heightfield Approximation: Shallow wave equation
Particle Systems: SPH

Chapter 15. Heightfield Approximations

默认读者都会

forall i,j do u[i,j] = u0[i,j]; v[i,j] = 0;
loop
forall i,j do v[i,j] += (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1])/4
- u[i,j];
forall i,j do v[i,j] *= 0.99;
forall i,j do u[i,j] += v[i,j];
endloop

Chapter 16. Smoothed Particle Hydrodynamics

smoothing kernel(常用poly6 kernel):

\[W_{\textrm{poly6}}(r)=\frac{315}{64\pi d^9}\left\{ \begin{aligned} &(d^2-r^2)^3&0\le r\le d\\ &0&\textrm{otherwise,} \end{aligned} \right. \]

compute a smooth density field:

\[\rho(x)=\sum_jm_jW(|x-x_j|) \]

Having the density values of the individual particles, we can now compute smoothed fields:

\[A_s(x)=\sum_jm_j\frac{A_j}{\rho_j}W(|x-x_j|) \]

\(A\): arbitrary attributes

\[\nabla A_s(x)=\sum_jm_j\frac{A_j}{\rho_j}\nabla W(|x-x_j|) \]

16.1 Pressure and Viscosity

更新密度：