这个系列是GAMES103：基于物理的计算机动画入门(GAMES 103: Intro to Physics Based Animation)的同步课程笔记。本课程是计算机动画的入门课程，着重介绍各种基于物理的动画仿真模拟技术。本节课主要介绍不可压缩流体的NS方程。
在上节课中我们简单介绍过Lagrangian方法和Eulerian方法的区别，简单来说Lagrangian方法是从粒子的角度进行分析，而Eulerian方法则是从流场所在空间的角度来分析。流体力学中的很多经典方法都是基于Eulerian方法来进行推导的，因此本节课会关注如何基于Eulerian方法进行流体仿真。
Grid Representation and Finite Differencing
Grid Representation
Eulerian方法的基本思想是对流场进行划分，其中最基本的方式是把整个流场划分成均匀分布的网格然后为每个网格分配相应的物理量(密度、速度、压力…)，这样流体的物理性质就可以通过离散后的场方程来描述。
Central Differencing
在网格的基础上我们还需要对场函数进行微分。参考上一节中Discretization部分，我们可以得到一维函数二阶精度的差分计算公式：
二维网格的计算方法与一维类似，只不过考虑每个格子上下左右四个相邻的格子进行差分。更进一步，我们可以把一阶差分定义在格子的边界上，进而计算每个格子中心处的二阶差分：
Discretized Laplacian
二阶差分的一个重要应用是用来构造离散Laplace算子(Laplacian operator)。对于连续函数，Laplace算子定义为函数在两个方向二阶偏导数的和。我们把这个定义推广到二维离散网格上即可获得离散Laplace算子：
Boundary Conditions
要求解微分方程还需要考虑边界条件。我们在上一节中已经介绍过两类边界条件：Dirichlet边界和Neumann边界。在二维网格上Dirichlet边界需要直接带入Laplace算子中，而Neumann边界则需要将边界点跳过：
Staggered Grid
在推导二阶差分时我们将一阶差分的函数值定义在格子的边界上，尽管这种做法仍然满足数学上的近似但却非常地不符合直觉。
为了解决这样的问题，我们可以把一些物理量定义在格子的边界而不是中心，这种做法称为staggered grid。在流体仿真时我们往往会把速度定义在格子的边界上：竖向的边界保存x轴方向的速度，横向的边界保存y方向的速度。这样的话每个格子中心处速度的导数就可以由格子边界的速度来进行推导。
除此之外，staggered grid另一个好处在于每个格子两个方向速度的变化量等于格子内流体体积的变化量。对于不可压缩流体我们需要保证在每个格子上的净流出量为0，即质量守恒(conservation of mass)。
Bilinear Interpolation
最后，如果需要计算流场中非中心位置处的物理量就需要通过插值来计算。对于二维网格，最常用的插值方法是双线性插值(bilinear interpolation)：我们需要提取\(\mathbf{x}\)附近的4个点并分别在两个方向上分别进行线性插值。而对于staggerred grid的情况，我们也需要提取四个相邻边界上的物理量进行插值。
Incompressible, Viscous Navier Stokes’ Equations
接下来我们介绍流体力学中最重要的方程，NS方程(Navier Stokes’ Equations)。NS方程本质是动量守恒方程，对于不可压缩流体NS方程可以化简为：
\[\frac{\partial \mathbf{u}}{\partial t} = \mathbf{g} - (\mathbf{u} \cdot \nabla) \mathbf{u} + \nu \Delta \mathbf{u} - \nabla p\]其中\(\mathbf{g}\)为外力引起的加速度，在大多数情况下等于重力加速度；\((\mathbf{u} \cdot \nabla) \mathbf{u}\)表示对流；\(\nu \Delta \mathbf{u}\)表示流体自身的扩散，\(\nu\)为流体的粘滞系数；\(\nabla p\)则是压强项。
除此之外流场还需要满足质量守恒：
\[\nabla \cdot \mathbf{u} = 0\]因此不可压缩流的仿真本质是求解这两个方程。显然要直接求解这两个微分方程是非常复杂的，在求解时我们可以考虑使用特征线法(method of characteristics)将整个方程拆分为若干个部分进行求解。对于NS方程，它可以拆分成以下4个部分：
- 利用外力\(\mathbf{g}\)更新速度
- 利用对流\((\mathbf{u} \cdot \nabla) \mathbf{u}\)更新速度
- 利用扩散\(\nu \Delta \mathbf{u}\)更新速度
- 利用压强\(\nabla p\)更新速度
External Acceleration
使用外力进行更新比较简单，只需要直接进行积分即可。而对于只有重力作为外力的情况则只需要更新y轴方向的速度即可。
Advection
关于对流项\((\mathbf{u} \cdot \nabla) \mathbf{u}\)我们可以把它进行展开得到加速度：
\[(\mathbf{u} \cdot \nabla) \mathbf{u} = \begin{bmatrix} u \cdot \frac{\partial u}{\partial x} + v \cdot \frac{\partial u}{\partial y} \\ u \cdot \frac{\partial v}{\partial x} + v \cdot \frac{\partial v}{\partial y} \\ \end{bmatrix}\]但如果直接进行积分的话则容易出现数值不稳定的情况。导致这种问题的原因在于计算对流项时的速度不应该是当前格子中的速度，从物理的角度上讲对流项中的速度是粒子运动到当前格子时携带的速度，因此进行积分时我们也应该使用粒子的速度而不是固定网格中的速度。
我们可以利用倒推的思路来解决这样的问题。对于格子边界上\(\mathbf{x}_0\)位置，它在上一时刻的位置可以近似表示为：
\[\mathbf{x}_1 = \mathbf{x}_0 - \Delta t \ \mathbf{u}(\mathbf{x}_0)\]然后通过插值来计算出上一时刻的速度\(\mathbf{u}(\mathbf{x}_1)\)来代替当前格子的速度\(\mathbf{u}(\mathbf{x}_0)\)即可。这样的方法称为semi-lagrangian method。
如果时间步长\(\Delta t\)比较大的话还可以采用逐步倒推的方式将\(\Delta t\)拆分成若干小步，然后依次向前递推粒子的位置并使用最后一步计算出的速度作为粒子在上一时刻的速度。
Diffusion
对于扩散项同样直接进行积分即可，但需要注意的是当\(\nu \Delta t\)比较大时同样会出现数值不稳定的问题。
要解决这个问题可以采用前面提过的缩减时间步长的方法，将\(\Delta t\)拆分成若干小步进行积分。
Pressure Projection
求解的最后一步是利用压强来对速度进行更新，对于网格来说同样非常简单。
那么如何计算网格上的压强呢？从物理本质上讲压强来自于流体的不可压缩性，因此我们需要使用质量守恒来更新压强：
整理后可以得到关于压强的方程组。我们只需要结合边界条件求解出压强再更更新速度场即可。
Air and Liquid
本节课最后简单介绍了烟尘和水面的模拟。对于烟尘而言，它的仿真过程可以大致分为2步：首先更新速度场，然后再利用semi-lagrangian来更新其它的物理量。
对于水的模拟则需要关注液体的不同表示方法。目前常用的表示方法包括volume-of-fluid和SDF两种，前者在进行仿真时会使用semi-lagrangian方法而后者则会使用level set方法。进行水的仿真时需要额外注意体积减少的问题，目前还没有完全解决的办法，一般需要额外对体积进行修正。
