流体模拟：基础

目录

• 流体模拟
  • 1.矢量微积分
    • 1.1梯度（Grant）
    • 1.2散度（Divergence）
    • 1.3旋度(Curl)
    • 1.4 拉普拉斯算子（Laplacian）
  • 2.Navier-Stokes偏微分方程
    • 2.1符号标记
    • 2.2动量方程（The Momentum Equation）
    • 2.3 拉格朗日描述与欧拉描述
    • 2.4不可压缩性
    • 2.5 丢弃粘度项
    • 2.6 边界条件
  • 3.N-S方程的分步求解
    • 3.1 时间步长
    • 3.2 网格结构
  • 4.对流算法
    • 4.1半拉格朗日对流算法（Semi-Lagrangian Advection）
    • 4.2 边界情况
    • 4.3时间步长大小
    • 4.4数值耗散
• 参考资料

流体模拟

1.矢量微积分

1.1梯度（Grant）

梯度实际上就是矢量的空间偏导数，且结果依然是一个矢量，2维的梯度如下

\[\nabla f(x, y)=\left(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}\right) \tag{1.1}\label{1.1} \]

有时也会采用如下形式来表示梯度：

\[\nabla f=\frac{\partial f}{\partial \vec{x}} \tag{1.2}\label{1.2} \]

有限微分形式(中心差分法)如下

\[\frac{f_{i+1, j}-f_{i-1, j}}{2 \delta x}, \frac{f_{i, j+1}-f_{i, j-1}}{2 \delta y} \tag{1.3}\label{1.3} \]

1.2散度（Divergence）

散度算子仅仅应用于向量场，它衡量在某一点处相应的矢量聚集或者发散程度，测量方向为径向，结果为标量。同时散度有重要的物理意义，它是“密度”离开给定空间区域的速率。

2维、3维形式的散度算子如下所示

\[\begin{gathered} \nabla \cdot \vec{u}=\nabla \cdot(u, v)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} \\ \nabla \cdot \vec{u}=\nabla \cdot(u, v, w)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z} \end{gathered} \tag{1.4}\label{1.4} \]

散度符号\(\nabla \cdot \vec{u}\)可以理解为梯度$\nabla \(与矢量\)\vec{u}$的点乘

\[\nabla \cdot \vec{u}=\left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}\right) \cdot(u, v, w)=\frac{\partial}{\partial x} u+\frac{\partial}{\partial y} v+\frac{\partial}{\partial z} w \tag{1.5}\label{1.5} \]

有限微分形式(中心差分法)：

\[\frac{u_{i+1, j}-u_{j-1, j}}{2 \delta x}+\frac{v_{i, j+1}-v_{i, j-1}}{2 \delta y} \tag{1.6}\label{1.6} \]

若矢量场散度为0，则称该矢量场无散度。在后面的Navier-Stokes方程中，他被应用于流速，他测量的是经过一小块流体表面的速率变化，让流体的的散度永远是零（即不可压缩状态）。

1.3旋度(Curl)

旋度衡量围绕某一点的旋转速度，测量方向为切向，三维形式的旋度是一个向量：

\[\nabla \times \vec{u}=\nabla \times(u, v, w)=\left(\frac{\partial w}{\partial y}-\frac{\partial v}{\partial z}, \frac{\partial u}{\partial z}-\frac{\partial w}{\partial x}, \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right) \tag{1.7}\label{1.7} \]

倒推到2维，我们取上式中的\(w=0\)，即矢量场为\((u,v,0)\)，2维向量场的旋度是一个标量:

\[\nabla \times \vec{u}=\nabla \times(u, v)=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} \tag{1.8}\label{1.8} \]

同样地,旋度符号\(\nabla \times \vec{u}\)我们可以理解为梯度\(\nabla\)与矢量场\(\vec{u}\)的叉乘:

\[\nabla \times \vec{u}=\left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}\right) \times(u, v, w) \tag{1.9}\label{1.9} \]

若矢量场旋度为0，则称该矢量场无旋度。若流体的旋度为0，即流体是不可旋转的状态。

1.4 拉普拉斯算子（Laplacian）

拉普拉斯算子定义为梯度的散度，符号表示为\(\nabla \cdot \nabla\),显然\(\nabla \cdot\)是散度,而后面的\(\nabla\)则为梯度，故拉普拉斯算子即梯度的散度，是一个二阶微分算子。2维、3维形式分别如下：

\[\begin{gathered} \nabla \cdot \nabla f=\frac{\partial^{2} f}{\partial x^{2}}+\frac{\partial^{2} f}{\partial y^{2}} \\ \nabla \cdot \nabla f=\frac{\partial^{2} f}{\partial x^{2}}+\frac{\partial^{2} f}{\partial y^{2}}+\frac{\partial^{2} f}{\partial z^{2}} \end{gathered} \tag{1.10}\label{1.10} \]

偏微分方程\(\nabla \cdot \nabla f=0\)被称为拉普拉斯方程;若右边为某个非0常数，即\(\nabla \cdot \nabla f=q\),如热传导方程，我们称之为泊松方程。更一般地，如果梯度再乘上一个标量a（如\(\frac {1}{\rho}\))，即\(\nabla \cdot (a\nabla f)=q\)，我们依旧称之为泊松问题。

有限微分形式(中心差分法)（2维形式）：

\[\frac{f_{i+1, i}-2 f_{i, j}+f_{i-1, j}}{(\delta x)^{2}}+\frac{f_{i, j+1}-2 f_{i, i}+f_{i, j-1}}{(\delta y)^{2}} \tag{1.11}\label{1.11} \]

如果网格单元是正方形（即\(\delta x=\delta y\)）,拉普拉斯算子可简化为：

\[\nabla^{2} p=\frac{f_{i+1, j}+f_{i-1, j}+f_{i, j+1}+f_{i, j-1}-4 f_{i, j}}{(\delta x)^{2}} \tag{1.12}\label{1.12} \]

2.Navier-Stokes偏微分方程

流体模拟器的构建主要围绕著名的不可压缩Navier−StokesNavier−Stokes方程展开，它是一个流体力学领域的偏微分方程，方程形式如下：

\[\begin{eqnarray*} \frac{\partial \vec{u}}{\partial t}+\vec{u} \cdot \nabla \vec{u}+\frac{1}{\rho} \nabla p=\vec{g}+\nu \nabla \cdot \nabla \vec{u} \tag{2.1} \label{2.1} \\ \nabla \cdot \vec{u}=0 \tag{2.2} \label{2.2} \end{eqnarray*} \]

现在我们将其剖析成一个个比较简单的部分来理解。

2.1符号标记

符号\(\vec {u}\)在流体力学中通常表示为流体的速度矢量，记3维的速度矢量\(\vec{u}=(u,v,w)\);

希腊字符\(\rho\)是流体的密度，对于水，该值大约为\(1000kg=m/s^3\)，而空气则大约为\(1.3kg=m/s^3\)

字符\(p\)代表压力，流体对任何物体施加的单位面积力

字符\(\vec{g}\)则是我们熟悉的重力加速度，通常取\((0,-9.81,0)m/s^3\)。我们约定y轴向上，而x轴和z轴在水平面上。另外补充一点，我们把其他的一些类似的力都累加\(\vec{g}\)上，也就是我们统一\(\vec{g}\)表示所有类似力之和，这类力我们称之为体积力（称之为体积力是因为它们的力是作用到整个流体而不只是流体的表面）；

希腊字符\(\vec{v}\)流体的运动粘度，它衡量流体的黏滞性。糖蜜之类的流体有非常高的粘度，而像酒精之类的流体有很低的粘度;

2.2动量方程（The Momentum Equation）

偏微分方程\(\eqref{2.1}\)我们称之为动量方程，它本质上就是我们熟悉的牛顿定律\(\vec{F}=m\vec{a}\)的形式，描述了施加在流体上的力是如何影响流体的运动。

假设我们用粒子系统来模拟流体，每个粒子代表流体的一小滴，每个粒子有各自的质量m、体积V和速度\(\vec{u}\)。为了让整个粒子系统运作起来，我们必须弄清楚每个粒子所承受的力的作用。牛顿第二定律告诉我们：\(\vec{F}=m\vec{a}\)，而根据加速度定义，我们有

\[\vec{a}=\frac{D\vec{u}}{Dt} \tag{2.3}\label{2.3} \]

符号\(D\)是指物质导数，所谓物质导数，就是对流体质点求导数，而且是针对流体质点（在这里就是流体粒子）而不是空间的固定点。因而牛顿第二定律就变成：

\[m\frac{D\vec{u}}{Dt}=\vec{F} \tag{2.4}\label{2.4} \]

流体粒子承受的力有：

• 一个最简单的力就是重力：\(m\vec{g}\)

• 其他的流体质点（或其他流体粒子）也会对当前的流体粒子产生作用力:

  流体内部的相互作用力之一便是压力，较大压力的区域会向较低压力区域产生作用力。值得注意的是，我们只关注施加在粒子上的压力的净合力。例如,若施加在粒子上压力在每个方向上都相等，那么它的压力的合力便为0。

  我们用压力的负梯度（取负是因为方向是由压力大的区域指向压力小的区域）来衡量在当前流体粒子处压力的不平衡性，即取\(-\nabla p\)。那么流体粒子所承受的压力就是对\(-\nabla p\)在整个流体粒子的体积上进行积分,为了简化，我们简单地将\(V\)与\(-\nabla p\)相乘，故粒子压力部分为\(-\nabla p\)。

• 其他的流体相互作用力则是由流体的黏性产生的：

  我们直观地把这种力理解为尽可能使得粒子以周围区域的平均速度移动的力，也就是使得粒子的速度与周围区域粒子速度的差距最小化

  拉普拉斯算子是衡量一个量与之周围区域该量平均值之差的算符，因而\(\nabla \cdot \nabla \vec{u}\)是当前粒子速度矢量与周围区域平均速度矢量之差。

  为了计算粘滞力，我们同样对\(\nabla \cdot \nabla \vec{u}\)整个粒子体积\(V\)上进行积分，与前面类似，我们简单取\(V\nabla \cdot \nabla \vec{u}\)。除此之外，我们还引进一个称为动力粘度系数的物理量，符号为\(\mu\)。因而粘滞力为\(V\mu\nabla \cdot \nabla \vec{u}\)

把重力、压力和粘滞力综合一起，我们可得：

\[m \frac{D \vec{u}}{D t}=\vec{F}=m \vec{g}-V \nabla p+V \mu \nabla \cdot \nabla \vec{u} \tag{2.5}\label{2.5} \]

当粒子系统中的粒子数量趋于无穷大，而每个粒子大小趋于0时，会产生一个问题：此时粒子的质量\(m\)和体积\(V\)变为0，此时上式子变得没有意义。为此，我们把\(\eqref{2.5}\)调整以下，两边同除以体积\(V\),又因\(\rho=m/V\),故有：

\[\rho \frac{D \vec{u}}{D t}=\rho \vec{g}- \nabla p+ \mu \nabla \cdot \nabla \vec{u} \tag{2.5}\label{2.6} \]

两边同时除以\(\rho\).同时定义运动粘度为\(v=\mu/\rho\),变成了

\[\frac{D \vec{u}}{D t}+ \frac{1}{\rho} \nabla p= \vec{g}+ v \nabla \cdot \nabla \vec{u} \tag{2.7}\label{2.7} \]

2.3 拉格朗日描述与欧拉描述

现在我们要把物质导数（Material derivative）\(\frac{D \vec{u}}{D t}\)弄清楚，我们需要了解两种描述方法：拉格朗日描述和欧拉描述。

当我们尝试研究流体或可变形固体的运动的时候，通常有两种方法来描述：拉格朗日描述（ Lagrangian viewpoint）、欧拉描述（Eulerian viewpoint）

• 拉格朗日描述（ Lagrangian viewpoint）

这种描述方法把物体看成是由类似于粒子系统的形式组成，固体或流体的每个点看作一个独立的粒子，粒子有各自相应的位置\(\vec{x}\)和速度\(\vec{u}\)。我们可以把粒子理解为组成物体的分子。对于我们通常采用拉格朗日描述法进行建模模拟，即用一系列离散的粒子集来构建，粒子之间通过网格相联系。

• 欧拉描述（Eulerian viewpoint）

  欧拉描述方法则采用了完全不同的角度，它常被用于流体力学中。与拉格朗日描述追踪每个物体粒子的方法不同，欧拉描述关注点是空间中的一个固定点，并考察在这个固定点上流体性质（如密度、速度、温度等）是如何随着时间变化的。流体流动经过这个固定点可能会导致这个固定点的物理性质发生一些变化（如一个温度较高的流体粒子流经这个固定点，后面紧跟着一个温度较低的流体粒子流过固定点，那么这个固定点的温度会降低，但是并没有任何一个流体粒子的温度发生了变化）。

用天气测量举个简单的例子：拉格朗日描述方法就是你乘坐在一个随风而飘的热气球上，测量周围空气的压力、密度和浑浊度等天气指标；而欧拉描述方法就是你固定在地面上，测量流过的空气的天气指标。

把拉格朗日描述法和欧拉描述法联系起来的关键点就是物质导数。

从拉格朗日描述法出发：假设有一群粒子，每个粒子都有各自的位置\(\vec{x}\)和速度\(\vec{v}\),记\(q\)为通用的物理量（如密度、速度和温度等），每个粒子有其对应的\(q\)值，方程\(q(t,\vec{x})\)描述在时间点\(t\)而位置\(\vec{x}\)的粒子对应的物理量值\(q\)。我们取对时间\(t\)的导数（注意用到了求导链式法则，以及\(\frac{\partial q}{\partial \vec{x}}=\nabla q\)和\(\vec{u}=\frac{d\vec{x}}{dt}\)):

\[\frac{d}{d t} q(t, \vec{x})=\frac{\partial q}{\partial t}+\nabla q \cdot \frac{d \vec{x}}{d t}=\frac{\partial q}{\partial t}+\nabla q \cdot \vec{u} \equiv \frac{D q}{D t} \tag{2.8}\label{2.8} \]

这就是物质导数。将其代入式子我们就得到了流体动量方程。物质导数针对的是流体质点在这里就是流体粒子）而不是空间的固定点，上式子写得完整些就是：

\[\frac{D q}{D t}=\frac{\partial q}{\partial t}+u \frac{\partial q}{\partial x}+v \frac{\partial q}{\partial y}+w \frac{\partial q}{\partial z} \tag{2.9}\label{2.9} \]

对于给定的速度场\(\vec{u}\),流体的物理性质如何在这个速度场\(\vec{u}\)下的变化计算我们称之为对流（advection）。一个最简单的对流方程，就是其物理量的物质导数为0，如下所示：

\[\frac{D q}{D t}=0 \Longrightarrow \frac{\partial q}{\partial t}+\vec{u} \cdot \nabla q=0 \tag{2.10}\label{2.10} \]

上式的意义即在拉格朗日视角观察下，每个流体粒子的物理量保持不变。

当物质导数作用于向量时，如：RGB,速度场时，一个简单方法就是分别对待每个分量。公式如下：

\[\frac{D \vec{u}}{D t}=\left[\begin{array}{c} D u / D t \\ D v / D t \\ D w / D t \end{array}\right]=\left[\begin{array}{c} \partial u / \partial t+\vec{u} \cdot \nabla u \\ \partial v / \partial t+\vec{u} \cdot \nabla v \\ \partial w / \partial t+\vec{u} \cdot \nabla w \end{array}\right]=\frac{\partial \vec{u}}{\partial t}+\vec{u} \cdot \nabla \vec{u} \tag{2.11}\label{2.11} \]

如果你想知道偏导数的具体细节：

\[\frac{D \vec{u}}{D t}=\left[\begin{array}{c} \frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z} \\ \frac{\partial v}{\partial t}+u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+w \frac{\partial v}{\partial z} \\ \frac{\partial w}{\partial t}+u \frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w \frac{\partial w}{\partial z} \end{array}\right] \tag{2.12}\label{2.12} \]

2.4不可压缩性

关于流体的压缩性的物理细节描述，可以看下3blue1brown的视频散度与旋度：麦克斯韦方程组、流体等所用到的语言

我们需要知道通常情况下流体的体积变化非常小（除开一些极端的情况，而且这些极端情况我们日常生活中较少出现）。可压缩流体的模拟涉及到非常复杂的情况，往往需要昂贵的计算资源开销，为此在计算机流体模拟中我们通常把所有的流体当作是不可压缩的，即它们的体积不会发生变化。

任取流体的一部分，设其体积为\(\Omega\)而其边界闭合曲面为\(\partial\Omega\),我们可以通过围绕边界曲面\(\partial\Omega\)对流体速度\(\vec{u}\)在曲面法线方向上的分量进行积分来衡量这块部分流体的体积变化速率：

\[\frac{d}{d t} \operatorname{Volume}(\Omega)=\iint_{\partial \Omega} \vec{u} \cdot n \tag{2.13}\label{2.13} \]

对于不可压缩的流体，其体积保持为某个常量，故其体积变化速率为0：

\[\iint_{\partial \Omega} \vec{u} \cdot n=0 \tag{2.14}\label{2.14} \]

由高斯散度定理（高斯公式），我们可以把上式转化为体积分：

\[\iint_{\partial \Omega} \vec{u} \cdot n=\iiint_{\Omega} \nabla \cdot \vec{u}=0 \tag{2.15}\label{2.15} \]

上式对任意的\(\Omega\)成立，即无论\(\Omega\)取何值，积分值均为0。这种情况下只有令积分函数取0方可成立，即对0积分无论\(\Omega\)取和值结果均为0。所以有：

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

这就是Navier−StokesNavier−Stokes方程中的不可压缩条件\(\eqref{2.2}\)。满足不可压缩条件的速度场被称为是无散度的，即在该速度场下流体体积既不膨胀也不坍缩，而是保持在一个常量。模拟不可压缩流体的关键部分就是使得流体的速度场保持无散度的状态，这也是流体内部压力的来源。

为了把压力与速度场的散度联系起来，我们在动量方程\(\eqref{2.1}\)两边同时取散度：

\[\nabla \cdot \frac{\partial \vec{u}}{\partial t}+\nabla \cdot(\vec{u} \cdot \nabla \vec{u})+\nabla \cdot \frac{1}{\rho} \nabla p=\nabla \cdot(\vec{g}+\nu \nabla \cdot \nabla \vec{u}) \tag{2.16}\label{2.16} \]

对于上式第一项，我们转变一下求导次序：

\[\frac{\partial}{\partial t}\nabla \cdot\vec{u} \tag{2.17}\label{2.17} \]

如果满足流体不可压缩条件，那么上式取值0（无散度），我们调整下上式可得关于压力的方程：

\[\nabla \cdot \frac{1}{\rho} \nabla p=\nabla \cdot(- \vec{u} \cdot \nabla \vec{u}+\vec{g}+\nu \nabla \cdot \nabla \vec{u}) \label{2.18} \tag{2.18} \]

2.5 丢弃粘度项

在某些流体如蜂蜜、小水珠等的模拟中，粘滞力起着非常重要的作用。但是在大多数流体动画模拟中，粘滞力的影响微乎其微，为此秉持着方程组越简单越好的原则，我们常常丢弃粘度项。当然这也不可避免地带来一些误差，事实上，在计算流体力学中尽可能地减少丢弃粘度项带来的误差是一个非常大的挑战。下面的叙述都是基于丢弃粘度项的前提。

丢弃了粘度项的Navier−StokesNavier−Stokes方程被称为欧拉方程，而这种理想的流体则是无粘度的。丢弃了粘度项的欧拉方程如下：

\[\begin{eqnarray*} \frac{D \vec{u}}{D t}+\frac{1}{\rho} \nabla p=\vec{g} \tag{2.19}\label{2.19} \\ \nabla \cdot \vec{u}=0 \tag{2.20}\label{2.20} \end{eqnarray*} \]

大多数的流体模拟的计算方程都是欧拉方程。

2.6 边界条件

目前为止我们讨论的都是流体内部的情况，然而边界部分也是流体模拟非常关键的部分。在流体模拟中我们仅仅关注两种边界条件：固体墙（solid walls）、自由面（free surfaces）。

固体墙顾名思义就是流体与固体接触的边界，用速度来描述很简单：流体既不会流进固体内部也不会从固体内部流出，因此流体在固体墙法线方向上的分量为0：

\[\vec{u} \cdot n=0 \tag{2.21}\label{2.21} \]

当然，上述是固体自身不移动的情况下。通常来说，流体速度在法线方向上的分量与固体的移动速度在法线方向上的分量应该保持一致:

\[\vec{v} \cdot n=\vec{u_{solid}} \cdot n \tag{2.22}\label{2.22} \]

上述的两个公式都是仅对流体速度在法线方向上的分量做了限制，对于无粘度的流体，切线方向上的流体速度与固体的移动速度无必然的联系。

自由面是另外一个非常重要的边界条件，它通常就是与另外一种流体相接壤的边界部分。例如在模拟水花四溅时，水流表面不与固体接触的都是自由面（如与空气这种流体接触）。因空气密度远小于水导致空气对水体的仿真影响非常小，为了简化模拟，我们将空气所占的空间设为某个固定大气压的区域，设为0是最方便的方案，此时自由面就是压强\(p=0\)的水体表面。

在小规模的流体仿真中，自由面的表面张力占据着非常重要的地位。在微观分子层面下，表面张力的存在是因为不同的分子相互吸引产生的力。从几何的角度来解释就是，表面张力就是促使流体的表面积尽可能小的一种力。物理学上，两种不同的流体之间实际上存在着与表面平均曲率成正比的压力骤变:

\[[p]=\lambda k \tag{2.23}\label{2.23} \]

上式\([p]\)记为压力之差。\(\lambda\)是表面张力系数，可以根据模拟的流体类型查找对应的张力系数（例如空气与水在室温下张力系数为\(\lambda \approx 0.073N/m\),。而\(k\)就是平均曲率，单位是\(m^{-1}\)。因为我们常常设空气的压力为0，因此水与空气交界的自由面的压力为：

\[p=\lambda k\tag{2.24}\label{2.24} \]

3.N-S方程的分步求解

了对以上对Navier−StokesNavier−Stokes方程的理论支撑，接下来我们就要如何用计算机来对该组偏微分方程进行离散化求解。为了程序的松耦合性以及使计算尽可能地高效、简单，在流体模型领域，我们将流体方程分成几个独立的步骤，然后按顺序先后推进。对于不可压缩的无粘度流体方程（即前面的欧拉方程），我们将其离散化成对流项（advection）\(\eqref{3.1}\)，体积力项（body force）\(\eqref{3.2}\)，压力/不可压缩项 \(\eqref{3.3}\)：

\[\begin{eqnarray*} \frac{D q}{D t}=0 \tag{3.1} \label{3.1} \\ \frac{\partial \vec{u}}{\partial t}=\vec{g} \tag{3.2} \label{3.2}\\ \left\{\begin{array}{l} \frac{\partial \vec{u}}{\partial t}+\frac{1}{\rho} \nabla p=0 \\ \nabla \cdot \vec{u}=0 \end{array}\right. \tag{3.3} \label{3.3} \end{eqnarray*} \]

• 对于对流项：我们用了一个通用量的符号\(q\)是因为我们不仅仅要对流体的速度进行对流，还需要对其他物理量进行对流。我们记对流项公式\(\eqref{3.1}\)的对流计算算法为\(advect(\vec{u},\triangle t,q)\)，即对于给定的时间步长\(\triangle t\)和速度\(\vec{u}\)，对物理量\(q\)进行对流。
• 对于体积力项：我们采用简单的前向欧拉法即可：\(\vec{u} \leftarrow \vec{u}+g \triangle t\)
• 对于压力/不可压缩项:我们用一个称为\(project(\triangle t,\vec{u})\)的算法，通过\(project(\triangle t,\vec{u})\)计算出正确的压力以确保速度场\(\vec{u}\)的无散度性质。欧拉方案不会着重研究具体粒子间的作用力，因而不会正向去求解\(\frac{1}{\rho} \nabla p\),它是利用流体不可压缩的特性，将速度场\(\vec{u}\)投影到散度为0的空间上，间接地解算了压力项。这种思想相当于，已知一个中间量\(\vec{u_{temp}}\),对这个中间量的唯一一个操作（如正向求解压力\(\frac{1}{\rho} \nabla p\)）不可行，但是知道最终量\(\vec{u_{final}}\)符号的一个性质（散度为0），于是只要将\(\vec{u_{temp}}\)投影到符合散度为0的特性平面上，即可间接地还原正向求解压力的操作，得到最终的速度场\(\vec{u_{temp}}\)。

对流项\(advect(\vec{u},\triangle t,q)\)的输入速度场\(\vec{u}\)要确保为无散度的状态，投影项\(project(\triangle t,\vec{u})\)确保流体体积保持不变，因而投影项输出的速度场必然是无散度。所以我们只要确保投影项\(project(\triangle t,\vec{u})\)输出的速度\(\vec{u}\)作为对流项\(advect(\vec{u},\triangle t,q)\)的输入即可，这时我们的分步求解流体方程的优势就体现出来了，其伪代码如下所示:

算法1：Fluid Simulation(\(\vec{u_n},\triangle t\))

1. 初始化速度场\(\vec{u_n}\),使得\(\vec{u_n}\)无散度
2. 对于每个时间步\(n=0,1,2...\)
3. 决定一个合理的时间步长\(\triangle t=t_{n+1}-t_n\)
4. 对流项计算$ \vec {u_A}= advect(\vec{u_n},\triangle t,q)$
5. 体积力项计算\(\vec {u_B}=project(\triangle t,\vec{u})\)
6. 无散度计算\(\vec {u_{n+1}}=project(\triangle t,\vec{u_b})\)

3.1 时间步长

在流体模拟算法中，确定适当的时间步长是算法的第一步。因为计算流体模拟的最后结果是呈现在屏幕上的，所以\(\triangle t\)的选取与屏幕的刷新率有重要的关系。若选取的\(\triangle t\)有\(t_n+\triangle t>t_{frame}\),那么必须做一个截断使\(\triangle t=t_{frame}-t_n\)。此外，流体模拟的三个步骤即对流项、体积力项、无散度投影项对时间步长\(\triangle t\)的要求不尽相同，要选择一个满足所有要求的最小时间步长能确保计算的收敛性。此外，一方面为了流体模拟的真实性，我们可能需要选取一个足够小的时间步长来复现流体的高质量细节。另一方面，有时高性能的需求又使得我们不能选取太小的时间步长去渲染一帧。假设一帧至少要进行三个时间步的模拟，那么\(\triangle t\)应该至少设成帧间隔时间的三分之一。

3.2 网格结构

欧拉法的整个流程都是基于网格的，所以合理的网格结构是算法高效的关键点。Harlow和Welch提出了一种经典的MAC（marker and cell）网格结构，许多不可压缩流体模拟的算法都在这个网格结构上呈现出了良好的效率。MAC网格是一种交叉排列的网格，不同类型的物理量被存储于网格的不同位置。以二维的网格为例，如下图所示，流体粒子的压力数据存储于网格的中心点\(P_{i,j}\)，而速度则沿着笛卡尔坐标被分成了两部分。水平方向的\(u\)成分被存储在了网格单元竖直边的中心处，例如网格单元\((i,j)\)和\((i+1,j)\)之间的水平速度记为\(u_{i+1/2,j}\)。垂直方向的\(v\)成分则被存储在了网格单元水平面的中心上。这样的存储方案十分有利于估算流体流进/流出某个网格单元的量。

扩展到三维的情况，MAC网格同样是交错排列的结构网格，如上图所示。压力数值存储在立方体网格单元的中心，三个速度分量分别被记录在立方体网格单元的三个表面的中心点上。在数值计算时，这样的分配方式使得我们可以准确地采用中心差分法计算压力梯度和速度的散度，同时克服了中心差分法的一个普遍的缺点。一维的情况为例，在网格顶点位置\(...,q_{i-1}，q_i,q_{i+1}...\)上估算量场\(q\)的导数，为了无偏（所谓无偏，就是不偏向左边或者右边）估计网格顶点i处的\(\frac{\partial q}{\partial x}\)，一种比较自然的方式就是采用一阶中心差分法:

\[\left(\frac{\partial q}{\partial x}\right)_{i} \approx \frac{q_{i+1}-q_{i-1}}{2 \Delta x} \tag{3.4} \label{3.4} \]

上式\(\eqref{3.4}\)是无偏的，且精确度为\(O(\triangle x^2)\)。而前向欧拉差分法偏向右边且精确度只有\(O(\triangle x)\)

\[\left(\frac{\partial q}{\partial x}\right)_{i} \approx \frac{q_{i+1}-q_{i}}{ \Delta x} \tag{3.5} \label{3.5} \]

然而，上式\(\eqref{3.4}\)存在着一个非常严重的问题：网格点\(i\)的估算导数完全忽略了\(q_i\)的值。数学上，只有常数函数的一阶导数为零。但是上式遇到了锯齿函数如\(q_i=(-1)^i\)时，它错误地将该类函数的导数估算为0，这种问题被称为零空间问题（null-space problem）。

交叉错排的MAC网格完美地克服了中心差分法的零空间问题，同时也保持了它的无偏二阶精度。在MAC网格上运用中心差分法，网格点i处的估算导数公式如下所示：

\[\left(\frac{\partial q}{\partial x}\right)_{i} \approx \frac{q_{i+1/2}-q_{i-1/2}}{ \Delta x} \tag{3.6} \label{3.6} \]

MAC网格确实给流体的压力计算和不可压缩性的处理带来了很大的便利，但与此同时也带来了一些其他方面的麻烦。如果我们要估算某个地方的速度向量，即便采样点恰好在网格点上我们也要做一些插值才能获取相应的速度向量。在网格点处，我们通常采用平均法，以二维为例：

\[\begin{gathered} \vec{u}_{i, j}=\left(\frac{u_{i-1 / 2, j}+u_{i+1 / 2, j}}{2}, \frac{v_{i, j-1 / 2}+v_{i, j+1 / 2}}{2}\right) \\ \vec{u}_{i+1 / 2, j}=\left(u_{i+1 / 2, j}, \frac{v_{i, j-1 / 2}+v_{i, j+1 / 2}+v_{i+1, j-1 / 2}+v_{i+1, j+1 / 2}}{4}\right) \\ \vec{u}_{i, j+1 / 2}=\left(\frac{u_{i-1 / 2, j}+u_{i+1 / 2, j}+u_{i-1 / 2, j+1}+u_{i+1 / 2, j+1}}{4}, v_{i, j+1 / 2}\right) \end{gathered} \tag{3.7} \label {3.7} \]

最后，在实现中下标索引一般没有浮点数之说，前面直接采用\(i+1/2\)的记法是为了便于叙述。一般约定如下：

\[\begin{gathered} p(i, j, k)=p_{i, j, k} \\ u(i, j, k)=u_{i-1 / 2, j, k} \\ v(i, j, k)=v_{i, j-1 / 2, k} \\ w(i, j, k)=w_{i, j, k-1 / 2} \end{gathered} \tag{3.8} \label{3.8} \]

而对于\(nx \times ny \times nz\)分辨率的网格，压力数值存储在\(nx \times ny \times nz\)的数组中，速度的\(u\)成分存储在\((nx+1) \times ny \times nz\)数组中，速度的\(v\)成分存储在\(nx \times (ny+1) \times nz\)数组中，速度的\(w\)成分存储在\(nx \times ny \times (nz+1)\)数组中。

4.对流算法

求解如下所示的对流方程是流体模拟的关键一步：

\[\frac{Dq}{Dt}=0 \tag{4.1} \label{4.1} \]

我们把这个对流数值计算的算法记为：

\[q^{n+1}=advect(\vec{u},\triangle t,q^n) \tag{4.2} \label{4.2} \]

公式\(\eqref{4.2}\)中的各个符号含义：

• \(\vec{u}\):在MAC网格上的离散化的速度场；
• \(\triangle t\):时间步长
• \(q^n\):当前的物理量场\(q\)（如流体密度、速度、燃烧物浓度等）；
• \(q^{n-1}\):经过对流后得到的新的量场。

在这里要特别注意，输入对流算法的速度\(\vec{u}\)必须是无散度的，否则模拟结果会出现一些奇怪的失真现象。

4.1半拉格朗日对流算法（Semi-Lagrangian Advection）

一维情况下，对流方程\(\eqref{4.1}\)写出偏微分的形式如下：

\[\frac{\partial q}{\partial t}+u\frac{\partial q} {\partial x}=0 \tag{4.3} \label{4.3} \]

分别采用前向欧拉差分法计算对时间的偏导和中心差分法计算对空间的偏导，我们有:

\[\frac{q_{i}^{n+1}-q_{i}^{n}}{\Delta t}+u_{i}^{n} \frac{q_{i+1}^{n}-q_{i-1}^{n}}{2 \Delta x}=0 \tag{4.4} \label {4.4} \]

转成以\(q^{n+1}_i\)维计算目标的显式公式，得：

\[q_{i}^{n+1}=q_{i}^{n}-\Delta t u_{i}^{n} \frac{q_{i+1}^{n}-q_{i-1}^{n}}{2 \Delta x} \tag{4.5} \label{4.5} \]

公式\(\eqref{4.1}\)有着非常严重的漏洞。

首先，前向欧拉法被证明是无条件不稳定的空间离散方法：无论取多么小Δ𝑡，随着时间步的推进，累积误差终将发散；然后标准中心差分方法不可避免地会出现的零空间问题，具有高频震荡性质的速度场对空间的导数被错误地计算为0或几乎为0，低离速度分量被分离出来，从而导致模拟效果中出现许多奇怪的高频摆动和震荡。

针对这些问题，研究者们提出了一个解然不同的、更加简单和更具物理直观意义的半拉格朗日法。之所以叫半拉格朗日法，是因为这种方法是以拉格朗日视角去解决欧拉视角的对流方程（“半”字的由来）

假设我们的目标室求解网格\(\vec{x_g}\)的在第n+1个时间步时关于物理量\(q\)的新值，记为\(q^{n+1}_G\)。在拉格朗日的视角下，我们可以寻找在第\(n+1\)时间步之前，是空间中的哪一个点上的流体粒子在速度场\(\vec{u}\)的作用下“流向”了\(\vec{x_G}\),我们记这个粒子在第\(n\)个时间步时的网格位置为\(\vec{x_p}\),则第\(n+1\)个时间步时\(q^{n+1}_G\)即为第\(n\)个时间步时\(\vec{x_p}\)的\(q^n_p\)。如下图为半拉格朗日对流法的示意图

半拉格朗日对流法的第一步就是要找\(\vec{x_p}\)，为此我们根据\(\vec{x_G}\)做反向的追踪。粒子位置对时间的导数就是速度场:

\[\frac{d\vec{x}}{d t}=\vec{u}(\vec{x}) \tag{4.6} \label{4.6} \]

经过一个时间步长\(\triangle t\)之后，粒子由\(\vec{x_p}\)移动到\(\vec{x_G}\)。为了得到\(\vec{x_p}\)，最简单的方法就是采用前向欧拉法进行倒推：

\[\vec{x_p}=\vec{x_G}-\triangle t\vec{u}（\vec{x_G}）\tag{4.7} \label{4.7} \]

然而前向欧拉法只有一阶的精度，若在不改变\(\triangle t\)的情况下提高精度，我们可以采用高阶的龙格库塔法（Runge-Kutta method）。采用二阶的龙格库塔法如下所示：

\[\begin{gathered} \vec{x}_{m i d}=\vec{x}_{G}-\frac{1}{2} \Delta t \vec{u}\left(\vec{x}_{G}\right) \\ \vec{x}_{P}=\vec{x}_{G}-\Delta t \vec{u}\left(\vec{x}_{m i d}\right) \end{gathered} \tag{4.7} \]

倒推得到\(\Delta t\)之前的网格位置\(\vec{x_P}\)一般不会恰好在网格顶点上，为此我们需要做些插值。三维模拟通常采用三线性插值，而二维的则采用双线性插值。

\[q_{G}^{n+1}=\text { interpolate }\left(q_{n}, \vec{x}_{P}\right) \tag{4.8} \label{4.8} \]

4.2 边界情况

若我们倒推得到\(\vec{x}_{P}\)仍然在流体的内部，那么做插值是完全没问题的。但\(\vec{x}_{P}\)在流体的边界之外呢？这种情况的出现的原因通常有两个：一个\(\vec{x}_{P}\)确确实实在流体的外部且即将流入流体内部，另一个是由前向欧拉法或龙格库塔法的数值计算方法带来的误差导致。

• 在一种情况下，我们应该知道当流体流入时其携带的物理量，此时我们将这个外部流入的物理量作为返回值即可。例如，第\(n\)个时间步时的外部流体以速度\(\vec{U}\)和温度\(T\)在第\(n+1\)个时间步时注入流体内\(\vec{x}_{G}\)的位置，那\(\vec{T}_G^{n+1}\)的值就为\(T\)。
• 在第二种由误差导致的情况下，一个适当的策略就是根据边界上的最近点外推出所求得物理量。在模拟某些流体时，外推变得很简单。例如，在模拟烟雾时我们简单地假设烟雾流体外部即空气的速度风场为某个常数\(\vec{U}\)（可能为0），这样边界上的速度场都取\(\vec{U}\)。但还有一些必须根据流体内部的已知量外推出未知量，这时情况就变得比较复杂了。具体如何外推将在后面介绍，目前我们只需要知道大概的步骤：首先寻找边界上的最近点，然后在最近点的领域内插值获取相应的物理量场。

4.3时间步长大小

对任何一种数值计算方法的主要的考虑点就是它是否稳定。幸运的是，半拉格朗日对流法已经被证明是一种无条件稳定的算法：无论\(\Delta t\)取多大，它永远不会出现数值爆炸的现象。因为每一个新值\(q\)的确定，都是通过对旧值得插值，无论是线性插值、双线性插值还是三线性插值，\(q\)的大小都是处于插值点之间，不会得到比原来插值点更大或者更小的值，因而\(q\)是有上下界的。这使得我们可以尽情地根据所需的模拟质量和模拟效率去调整时间步长。

但是在实践中，时间步长的大小也不能选得太过极端，否则会产生一些奇观的现象。Foster和Fekiw提出了一个对\(\Delta t\)的限制：流体粒子在ΔtΔt内的倒推轨迹最多经过某个常数个网格单元为宜，例如5个：

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

公式\(\eqref{4.9}\)中，\(u_{max}\)是速度场的最大值，我们可以简单地取 存储在网格中的最大速度值。一个更鲁棒的方法考虑了体积力（如重力、浮力等）对最大速度的影响：

\[u_{\max }=\max \left(\left|u^{n}\right|\right)+\Delta t|g| \tag{4.10} \label {4.10} \]

将不等式(4.9)(4.9)的最大值带入公式\(\eqref{4.10}\)，我们有：

\[u_{\max }=\max \left(\left|u^{n}\right|\right)+\frac{5 \Delta x}{u_{\max }}|g| \tag{4.11} \label {4.11} \]

取一个简单的速度上界（简化了公式\(\eqref{4.11}\)），\(u_{max}\):

\[u_{\max }=\max \left(\left|u^{n}\right|\right)+\sqrt{5 \Delta x g} \tag{4.12} \label {4.12} \]

这样确保了\(u_{max}\)始终为正，且避免公式\(\eqref{4.9}\)的除0错误。

关于时间步长的讨论离不开CFL（以Courant、Friedrichs、Lewy三人的名字命名）条件。CFL条件是一个简单而直观的判断计算是否收敛的必要条件。它的直观物理解释就是时间推进求解的速度必须大于物理扰动传播的速度，只有这样才能将物理上所有的扰动俘获到。满足CFL条件意味着当ΔxΔx和ΔtΔt趋于取极限00时，数值计算所求的解就会收敛到原微分方程的解。

对于半拉格朗日对流法，其满足CFL条件当且仅当在极限情况下，追踪得到的粒子轨迹足够逼近真实的轨迹。足够逼近的意思是经过正确的网格插值能够得到正确的依赖域（即差分格式的依赖域包含了原微分方程的依赖域），追踪的轨迹就会收敛到正确真实的轨迹。

因而，对于采用标准的显式有限差分法的对流方程求解，为了保证收敛，我们要求\(q^{n+1}\)的新值是由以当前网格点为中心、以\(C\Delta x\)（\(C\)是一个小的整数常量）为半径的邻域范围内插值得到：

\[\Delta t \leq C \frac{\Delta x}{|\vec{u}|} \tag{4.13} \label {4.13} \]

公式\(\eqref{4.13}\)中的\(C\)被称为CFL数，因而不等式(4.9)(4.9)可以看成是公式\(\eqref{4.13}\)取CFL数为5得到。

4.4数值耗散

对流算法在对流获取新的物理量场\(q^{n+1}_i\)时会进行一些插值操作，插值不可避免地会平滑物理量场，这带来了一些数值耗散。一次两次的数值耗散不会由太大的影响，但是在流体模拟中我们会在每个时间步都进行对流运算，反反复复的平滑操作将数值耗散不断扩大，损失大量的流体细节。

以一维的对流项计算为例，流体速度为常量\(U>0\):

\[\frac{\partial q}{\partial t}+u \frac{\partial q}{\partial x}=0 \tag{4.14} \label {4.14} \]

假设\(\Delta t < \frac{\Delta x}{u}\),即单个时间步长内粒子追踪轨迹小于单个网格单元的大小。我们的目标点时\(x_i\),则倒推地道道阿粒子位置就落在了\({x_{i-1},x_i}\)上的\(x_i-\Delta t u\),r然后进行线性插值得到\(q^{n+1}_i\):

\[q^{n+1}=\frac{\Delta t u}{\Delta x} q_{i-1}^{n}+\left(1-\frac{\Delta t u}{\Delta x}\right) q_{i}^{n} \tag{4.15} \label {4.15} \]

将公式\(\eqref{4.15}\)整理一下，有：

\[q_{i}^{n+1}=q_{i}^{n}-\Delta t u \frac{q_{i}^{n}-q_{i-1}^{n}}{\Delta x} \tag{4.16} \label {4.16} \]

公式\(\eqref{4.16}\)实际上正好就是采用时间上的前向欧拉差分法和空间上的单向有限差分法的欧拉方案，把\(q^n_i\)看成是\(q^n\)关于\(x_i\)的函数，对\(q^n_{i-1}\)进行泰勒级数展开：

\[q_{i-1}^{n}=q_{i}^{n}-\left(\frac{\partial q}{\partial x}\right)_{i}^{n} \Delta x+\left(\frac{\partial^{2} q}{\partial x^{2}}\right)_{i}^{n} \frac{\Delta x^{2}}{2}+O\left(\Delta x^{3}\right) \tag{4.17} \label {4.17} \]

将公式\(\eqref{4.17}\)代入公式\(\eqref{4.16}\)，并做一些变量消去，可得：

\[q_{i}^{n+1}=q_{i}^{n}-\Delta t u\left(\frac{\partial q}{\partial x}\right)_{i}^{n}+\Delta t u \Delta x\left(\frac{\partial^{2} q}{\partial x^{2}}\right)_{i}^{n}+O\left(\Delta x^{2}\right) \tag{4.18} \label {4.18} \]

在二阶截断误差的情况下，结合公式\(\eqref{4.18}\)和公式\(\eqref{4.14}\)，有：

\[\frac{\partial q}{\partial t}+u \frac{\partial q}{\partial x}=u \Delta x\left(\frac{\partial^{2} q}{\partial x^{2}}\right) \tag{4.19} \label {4.19} \]

右边就是对流方程计算时引入的额外类似粘度乘上系数\(u\Delta x\)的项。这也就是说，当我们采用简单的半拉格朗日法去求解无粘度的对流方程时，模拟的结果却看起来我们像时在模拟有粘度的流体。这就是数值耗散！当然，当\(\Delta x \rightarrow 0\)时，这个数值耗散系数也会趋于00，所以取时间步无穷小时能够得到正确的模拟结果，但这需要耗费巨额的计算资源开销。我们通常模拟的流体大多数都是无粘度的，所以如何减少这个数值耗散是个至关重要的难题。

一个简单有效的修复数值耗散的方法就是采用更加锐利的插值方法，从而尽可能地减少由插值带来的数值耗散。在一维的情况时，我们采用三次插值（cubic interpolant）如下公式\(\eqref{4.21}\)，而不是简单的一次线性插值\(\eqref{4.20}\)：

\[\begin{eqnarray*} q \approx(1-s) x_{i}+s x_{i+1} \tag{4.20} \label {4.20}\\ q \approx\left[-\frac{1}{3} s+\frac{1}{2} s^{2}-\frac{1}{6} s^{3}\right] q_{i-1}+\left[1-s^{2}+\frac{1}{2}\left(s^{3}-s\right)\right] q_{i} \tag{4.21} \label {4.21}\\ +\left[s+\frac{1}{2}\left(s^{2}-s^{3}\right)\right] q_{i+1}+\left[\frac{1}{6}\left(s^{3}-s\right)\right] q_{i+2} \end{eqnarray*} \]

扩展到二维或者三维就是双三次插值（bicubic interpolation）或三三次插值（tricubic interpolation）。以二维情况为例，我们可以先沿着\(x\)轴做第一遍的三次插值如公式\(\eqref{4.22}\)，然后再沿着\(y\)轴做第二遍插值如公式\(\eqref{4.23}\)：

\[\begin{eqnarray*} q_{j-1}=w_{-1}(s) q_{i-1, j-1}+w_{0}(s)+q_{i, j-1}+w_{1}(s) q_{i+1, j-1}+w_{2}(s) q_{i+2, j-1}, \tag{4.22} \label {4.22}\\ q_{j}=w_{-1}(s) q_{i-1, j}+w_{0}(s)+q_{i, j}+w_{1}(s) q_{i+1, j}+w_{2}(s) q_{i+2, j} \\ q_{j+1}=w_{-1}(s) q_{i-1, j+1}+w_{0}(s)+q_{i, j+1}+w_{1}(s) q_{i+1, j+1}+w_{2}(s) q_{i+2, j+1} \\ q_{j+2}=w_{-1}(s) q_{i-1, j+2}+w_{0}(s)+q_{i, j+2}+w_{1}(s) q_{i+1, j+2}+w_{2}(s) q_{i+2, j+2} \\ q=w_{-1}(t) q_{j-1}+w_{0}(t) q_{j}+w_{1}(t) q_{j+1}+w_{2}(t) q_{j+2} \tag{4.23} \label {4.23} \end{eqnarray*} \]

当然也可以先沿着\(y\)轴，然后再沿着\(x\)轴做插值操作。

参考资料

有限差分法（FDM）：微分方程数值求解——有限差分法

文献：《FLUID SIMULATION SIGGRAPH 2007 Course Notes》,《GPU Games》

博客：流体模拟Fluid Simulation：流体模拟基础

视频：散度与旋度：麦克斯韦方程组、流体等所用到的语言