GAMES103课程笔记11-Eulerian Fluids

这个系列是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}\]

但如果直接进行积分的话则容易出现数值不稳定的情况。导致这种问题的原因在于计算对流项时的速度不应该是当前格子中的速度，从物理的角度上讲对流项中的速度是粒子运动到当前格子时携带的速度，因此进行积分时我们也应该使用粒子的速度而不是固定网格中的速度。