流体模拟原理与实现 本文来源：YangWC's Blog 《流体模拟Fluid Simulation：流体模拟基础》本文主要参考文献《FLUID SIMULATION SIGGRAPH 2007 Course Notes》。一、矢量微积分  高等数学中太多数讨论的是一维的微积分，而矢量微积分则是一维微积分的高维扩展。矢量微积分的三个基础算子：梯度（符号为 $ \nabla $ ），散度（符号为 $ \nabla \cdot $ )，旋度（符号为 $ \nabla \times $ ），在此基础上流体力学中经常用到的还有拉普拉斯算子。1、梯度（Gradient）  梯度实际上就是矢量的空间偏导数，且结果依然是一个矢量，$2$维的梯度如下：$$ ∇f(x,y)=(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}) \tag {1.1}$$  依此类推，$3$维的梯度有如下形式：$$∇f(x,y,z)=(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z}) \tag {1.2}$$  有时也会采用如下形式来表示梯度：$$∇f=\frac{\partial f}{\partial \vec x} \tag {1.3}$$  梯度通常用来近似计算函数值（实际上就是一维形式的推广)：$$f(\vec x+\Delta \vec x)\approx f(\vec x)+∇f(\vec x)\cdot \Delta \vec x \tag {1.4}$$  同样的，多个函数的梯度就构成了一个矩阵：$$∇\vec F=∇(f,g,h)=\left( \begin{matrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} & \frac{\partial g}{\partial z} \\ \frac{\partial h}{\partial x} & \frac{\partial h}{\partial y} & \frac{\partial h}{\partial z} \\ \end{matrix} \right) =\left( \begin{matrix}∇f\\ ∇g\\ ∇h\\ \end{matrix} \right) \tag {1.5}$$2、散度（Divergence）  散度算子仅仅应用于向量场，它衡量在某一点出相应的矢量聚集或者发散程度，测量方向为径向，结果为标量。$2$维、$3$维形式的散度算子如下所示：$$∇\cdot \vec u=∇\cdot (u,v)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}$$$$∇\cdot \vec u=∇\cdot (u,v,w)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z} \tag {1.6}$$  输入是矢量，而输出为标量。类比梯度，散度符号$∇\cdot \vec u$可以理解为梯度$∇$与矢量$\vec u$的点乘：$$∇\cdot \vec u=(\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z})\cdot (u,v,w)=\frac{\partial}{\partial x}u+\frac{\partial}{\partial y}v+\frac{\partial}{\partial z}w \tag {1.7}$$  若矢量场散度为$0$，则称该矢量场无散度。3、旋度（Curl）  旋度衡量围绕某一点的旋转速度，测量方向为切向，三维形式的旋度是一个向量：$$∇\times \vec u=∇\times (u,v,w) =(\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}) \tag {1.8}$$  倒推到2维，我们取上式中的$w=0$，即矢量场为$(u,v,0)$，2维向量场的旋度是一个标量：$$∇\times \vec u=∇\times (u,v)=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} \tag {1.9}$$  同样地，旋度符号$∇\times \vec u$我们可以理解为梯度$∇$与矢量$\vec u$的叉乘：$$∇\times \vec u= (\frac{\partial }{\partial x}, \frac{\partial }{\partial y}, \frac{\partial }{\partial z})\times(u,v,w) \tag {1.10}$$  若矢量场旋度为0，则称该矢量场无旋度。4、拉普拉斯算子（Laplacian）  拉普拉斯算子定义为梯度的散度，符号表示为$∇\cdot∇$，显然$∇\cdot$是散度，而后面的$∇$则为梯度，故拉普拉斯算子即梯度的散度，是一个二阶微分算子。2维、3维形式分别如下：$$∇\cdot∇f=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2}$$$$∇\cdot∇f=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2}+\frac{\partial^2f}{\partial z^2} \tag {1.11}$$  简言之，拉普拉斯算子定义如下：$$∇\cdot∇f=\Sigma_{i=1}^n\frac{\partial^2f}{\partial x_i^2} \tag {1.12}$$  偏微分方程$∇\cdot ∇f=0$被称为拉普拉斯方程；而如果右边为某个非0常数，即$∇\cdot ∇f=q$，我们称之为泊松方程。更一般地，如果梯度再乘上一个标量$a$（如$1/ρ$)，即$∇\cdot (a∇f)=q$，我们依旧称之为泊松问题。二、Naiver−StokesNaiver−Stokes偏微分方程组  流体模拟器的构建主要围绕著名的不可压缩Navier−StokesNavier−Stokes方程展开，它是一个流体力学领域的偏微分方程，方程形式如下：$$\frac{\partial \vec u}{\partial t}+\vec u\cdot ∇\vec u+\frac1\rho∇p=\vec g+\nu∇\cdot∇\vec u \tag {2.1}$$$$∇\cdot\vec u=0 \tag {2.2}$$1、符号标记  我们有必要定义一些物理量的符号用以标记：  符号$\vec u$在流体力学中通常表示为流体的速度矢量，记$3$维的速度矢量$\vec u=(u,v,w)$；  希腊字符$\rho$是流体的密度，对于水，该值大约为$1000kg/m^3$，而空气则大约为$1.3kg/m^3$；  字符$p$代表压力，流体对任何物体施加的单位面积力；  字符$\vec g$则是我们熟悉的重力加速度，通常取$(0,-9.81,0)m/s^2$。我们约定$y$轴向上，而$x$轴和$z$轴在水平面上。另外补充一点，我们把其他的一些类似的力都累加到$\vec g$上，也就是我们统一用$\vec g$表示所有类似力之和，这类力我们称之为体积力（称之为体积力是因为它们的力是作用到整个流体而不只是流体的表面）；  希腊字符$\nu$是流体的运动粘度，它衡量流体的黏滞性。糖蜜之类的流体有非常高的粘度，而像酒精之类的流体有很低的粘度；  其它一些矢量微积分的符号算子前面已经提到过，不再赘述。2、动量方程  偏微分方程$(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}$$  符号$D$是指物质导数，所谓物质导数，就是对流体质点求导数，而且是针对流体质点（在这里就是流体粒子）而不是空间的固定点。因而牛顿第二定律就变成：$$m\frac{D\vec u}{Dt}=\vec F \tag {2.4}$$  那么流体粒子承受的力有哪些呢？一个最简单的力就是重力：$m\vec g$。而其他的流体质点（或其他流体粒子）也会对当前的流体粒子产生作用力。流体内部的相互作用力之一便是压力，较大压力的区域会向较低压力区域产生作用力。值得注意的是，我们只关注施加在粒子上的压力的净合力。例如，若施加在粒子上压力在每个方向上都相等，那么它的压力的合力便为0。我们用压力的负梯度（取负是因为方向是由压力大的区域指向压力小的区域）来衡量在当前流体粒子处压力的不平衡性，即取$-∇p$。那么流体粒子所承受的压力就是对$-∇p$在整个流体粒子的体积上进行积分，为了简化，我们简单地将$V$与$-∇p$相乘，故粒子压力部分为$-V∇p$。  其他的流体相互作用力则是由流体的黏性产生的，我们直观地把这种力理解为尽可能使得粒子以周围区域的平均速度移动的力，也就是使得粒子的速度与周围区域粒子速度的差距最小化。拉普拉斯算子是衡量一个量与之周围区域该量平均值之差的算符，因而$∇\cdot∇\vec u$是当前粒子速度矢量与周围区域平均速度矢量之差。为了计算粘滞力，我们同样对$∇\cdot∇\vec u$在整个粒子体积$V$上进行积分，与前面类似，我们简单取$V∇\cdot∇\vec u$。除此之外，我们还引进一个称为动力粘度系数的物理量，符号为$\mu$。因而粘滞力为$V\mu∇\cdot∇\vec u$。  把重力、压力和粘滞力综合一起，我们可得：$$m\frac{D\vec u}{Dt}=\vec F=m\vec g-V∇p+V\mu∇\cdot∇\vec u \tag {2.5}$$  当粒子系统中的粒子数量趋于无穷大，而每个粒子大小趋于$0$时，会产生一个问题：此时每个粒子的质量$m$和体积$V$变为$0$，此时上式变得没有意义。为此，我们把$(2.5)$式调整一下，两边同除以体积$V$，又因$\rho=m/V$，故有：$$\rho\frac{D\vec u}{Dt}=\rho\vec g-∇p+\mu∇\cdot∇\vec u \tag {2.6}$$  两边同除以$\rho$，移项调整：$$\frac{D\vec u}{Dt}+\frac1\rho∇p=\vec g+\frac\mu\rho∇\cdot∇\vec u \tag {2.7}$$  为了进一步简化，定义运动粘度为$\nu=\mu/\rho$，式$(2.7)$变为：$$\frac{D\vec u}{Dt}+\frac1\rho∇p=\vec g+\nu∇\cdot∇\vec u \tag {2.8}$$  我们已经快把动量方程推导出来，现在我们要把物质导数$\frac{D\vec u}{Dt}$弄清楚，为此，我们需要了解两种描述方法：拉格朗日描述和欧拉描述。3、拉格朗日描述与欧拉描述  当我们尝试研究流体或可变形固体的运动的时候，通常有两种方法来描述：拉格朗日描述（ Lagrangian viewpoint）、欧拉描述（Eulerian viewpoint）。  拉格朗日描述方法是我们比较熟悉的方法，这种描述方法把物体看成是由类似于粒子系统的形式组成，固体或流体的每个点看作一个独立的粒子，粒子有各自相应的位置$\vec x$和速度$\vec u$。我们可以把粒子理解为组成物体的分子。对于我们通常采用拉格朗日描述法进行建模模拟，即用一系列离散的粒子集来构建，粒子之间通过网格相联系。  欧拉描述方法则采用了完全不同的角度，它常被用于流体力学中。与拉格朗日描述追踪每个物体粒子的方法不同，欧拉描述关注点是空间中的一个固定点，并考察在这个固定点上流体性质（如密度、速度、温度等）是如何随着时间变化的。流体流动经过这个固定点可能会导致这个固定点的物理性质发生一些变化（如一个温度较高的流体粒子流经这个固定点，后面紧跟着一个温度较低的流体粒子流过固定点，那么这个固定点的温度会降低，但是并没有任何一个流体粒子的温度发生了变化）。  用天气测量举个简单的例子：拉格朗日描述方法就是你乘坐在一个随风而飘的热气球上，测量周围空气的压力、密度和浑浊度等天气指标；而欧拉描述方法就是你固定在地面上，测量流过的空气的天气指标。  欧拉描述法似乎看起来带来了一些不必要的复杂度，但是目前大多数的流体模拟器都是基于欧拉描述法，这是因为欧拉描述法相比于拉格朗日描述法有一些不可比拟的优点：欧拉描述法能够更加方便地计算一些物理量的空间导数（例如压力梯度和粘度）；而如果用粒子方法的话（即拉格朗日描述法），那么计算物理量相对于空间位置的变化是比较难的。  把拉格朗日描述法和欧拉描述法联系起来的关键点就是物质导数。首先从拉格朗日描述法出发，假设有一群粒子，每个粒子都有各自的位置$\vec x$和速度$\vec u$。记$q$为通用的物理量（如密度、速度和温度等），每个粒子有其对应的$q$值。方程$q(t,\vec x)$描述在时间点$t$而位置为$\vec x$的粒子对应的物理量值$q$。则一个粒子的物理量$q$随时间$t$的变化率是多少？这是一个拉格朗日描述角度下的问题，我们取对时间$t$的导数（注意用到了求导链式法则，以及$\frac{\partial q}{\partial \vec x}=∇q$和$\vec u=\frac{d\vec x}{dt}）$：$$\frac d{dt}q(t,\vec x)=\frac{\partial q}{\partial t}+∇q\cdot\frac{d\vec x}{dt}=\frac{\partial q}{\partial t}+∇q\cdot\vec u\equiv\frac{Dq}{Dt} \tag {2.9}$$  这就是物质导数。把式$(2.9)$代入式$(2.8)$我们就得到了流体动量方程$(2.1)$。物质导数针对的是流体质点（在这里就是流体粒子）而不是空间的固定点。式$(2.9)$写完整一点就是：$$ \frac{Dq}{Dt}=\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.10} $$  对于给定的速度场$\vec u$， 流体的物理性质如何在这个速度场$\vec u$下变化的计算我们称之为对流（advection）。一个最简单的对流方程，就是其物理量的物质导数为$0$，如下所示：$$\frac{Dq}{Dt}=0\implies\frac{\partial q}{\partial t}+\vec u\cdot ∇q = 0 \tag {2.11}$$  公式$(2.11)$的意义即在拉格朗日视角观察下，每个流体粒子的物理量保持不变。4、不可压缩性  关于流体的压缩性在此不做过多的物理细节描述，只需知道一点：通常情况下流体的体积变化非常小（除开一些极端的情况，而且这些极端情况我们日常生活中较少出现）。可压缩流体的模拟涉及到非常复杂的情况，往往需要昂贵的计算资源开销，为此在计算机流体模拟中我们通常把所有的流体当作是不可压缩的，即它们的体积不会发生变化。  任取流体的一部分，设其体积为$\Omega$而其边界闭合曲面为$\partial\Omega$，我们可以通过围绕边界曲面$\partial\Omega$对流体速度$\vec u$在曲面法线方向上的分量进行积分来衡量这块部分流体的体积变化速率：$$\frac d{dt}Volume(\Omega)=\int\int_{\partial\Omega}\vec u\cdot n \tag{2.12}$$  对于不可压缩的流体，其体积保持为某个常量，故其体积变化速率为$0$：$$\int\int_{\partial\Omega}\vec u\cdot n=0 \tag {2.13}$$  由高斯散度定理，我们可以把式$(2.13)$转换为体积分：$$\int\int_{\partial\Omega}\vec u\cdot n=\int\int\int_\Omega∇\cdot \vec u=0 \tag{2.14}$$  式$(13)$应该对任意的$\Omega$成立，意即无论$\Omega$取何值，积分值均为$0$。这种情况下只有令积分函数值取$0$方可成立，即对$0$积分无论$\Omega$取何值结果均为$0$。所以有：$$∇\cdot \vec u=0 \tag{2.15}$$  这就是$Navier-Stokes$方程中的不可压缩条件$(2.2)$。满足不可压缩条件的速度场被称为是无散度的，即在该速度场下流体体积既不膨胀也不坍缩，而是保持在一个常量。模拟不可压缩流体的关键部分就是使得流体的速度场保持无散度的状态，这也是流体内部压力的来源。  为了把压力与速度场的散度联系起来，我们在动量方程$(2.1)$两边同时取散度：$$ ∇\cdot\frac{\partial \vec u}{\partial t}+∇\cdot(\vec u\cdot ∇\vec u)+∇\cdot\frac1\rho∇p=∇\cdot(\vec g+\nu∇\cdot∇\vec u) \tag {2.16} $$  对于上式$(2.16)$第一项，我们转变一下求导次序：$$\frac {\partial}{\partial t}∇\cdot\vec u \tag {2.17}$$  如果满足流体不可压缩条件，那么式$(2.17)$取值$0$（因为无散度），然后我们调整一下式$(2.16)$可得关于压力的方程：$$∇\cdot\frac1\rho∇p=∇\cdot(-\vec u\cdot ∇\vec u+\vec g+\nu∇\cdot∇\vec u) \tag{2.18}$$5、丢弃粘度项  在某些流体如蜂蜜、小水珠等的模拟中，粘滞力起着非常重要的作用。但是在大多数流体动画模拟中，粘滞力的影响微乎其微，为此秉持着方程组越简单越好的原则，我们常常丢弃粘度项。当然这也不可避免地带来一些误差，事实上，在计算流体力学中尽可能地减少丢弃粘度项带来的误差是一个非常大的挑战。下面的叙述都是基于丢弃粘度项的前提。  丢弃了粘度项的$Navier-Stokes$方程被称为欧拉方程，而这种理想的流体则是无粘度的。丢弃了粘度项的欧拉方程如下：$$\frac{D\vec u}{Dt}+\frac1\rho∇p=\vec g \tag {2.19}$$$$∇\cdot\vec u=0 \tag{2.20}$$  大多数的流体模拟的计算方程都是欧拉方程。6、边界条件  目前为止我们讨论的都是流体内部的情况，然而边界部分也是流体模拟非常关键的部分。在流体模拟中我们仅仅关注两种边界条件：固体墙（solid walls）、自由面（free surfaces）。  固体墙顾名思义就是流体与固体接触的边界，用速度来描述很简单：流体既不会流进固体内部也不会从固体内部流出，因此流体在固体墙法线方向上的分量为$0$：$$\vec u\cdot n=0 \tag {2.21}$$  当然，上述是固体自身不移动的情况下。通常来说，流体速度在法线方向上的分量与固体的移动速度在法线方向上的分量应该保持一致：$$\vec u\cdot n=\vec u_{solid}\cdot n \tag{2.22}$$  上述的两个公式都是仅对流体速度在法线方向上的分量做了限制，对于无粘度的流体，切线方向上的流体速度与固体的移动速度无必然的联系。  自由面是另外一个非常重要的边界条件，它通常就是与另外一种流体相接壤的边界部分。例如在模拟水花四溅时，水流表面不与固体接触的都是自由面（如与空气这种流体接触）。因空气密度远小于水导致空气对水体的仿真影响非常小，为了简化模拟，我们将空气所占的空间设为某个固定大气压的区域，设为$0$是最方便的方案，此时自由面就是压强$p=0$的水体表面。  在小规模的流体仿真中，自由面的表面张力占据着非常重要的地位。在微观分子层面下，表面张力的存在是因为不同的分子相互吸引产生的力。从几何的角度来解释就是，表面张力就是促使流体的表面积尽可能小的一种力。物理学上，两种不同的流体之间实际上存在着与表面平均曲率成正比的压力骤变：$$ [p]=\lambda k. \tag {2.23} $$  公式$(2.23)$中的$[p]$记为压力之差。$\lambda$是表面张力系数，可以根据模拟的流体类型查找对应的张力系数（例如空气与水在室温下张力系数为$\lambda \approx 0.073N/m$）。而$k$就是平均曲率，单位为$m^{-1}$。又因为我们常常设空气的压力为$0$，因此水与空气交界的自由面的压力为：$$p=\lambda k \tag {2.24}$$三、N-S方程的分步求解1、数值近似  上面已经基本很详细的剖析了NS方程各项的意义和背后的数学建模原理，以及针对常见流体，也就是非粘性不可压缩流体的简化形式， 然而即使经过简化，变成了欧拉方程，解决起来依然不是很简单，无法求出解析解，因此流体仿真一般采用数值近似的方式。  数值近似直接解决整个欧拉方程是很复杂的，因此一般采取分治法的思想，分而治之，把欧拉方程分离成几个部分分别解决， 然后再结合起来，这在数学上是合理的，假设：$$\frac{dq}{dt}=f(q) + g(q)$$  这里f和g函数可以是任意的黑盒子函数代表分离的部分，我们可以将这个微分方程使用显式欧拉分离成：$$\tilde{q} =q^{n} + \Delta t f(q^{n})$$$$q^{n+1} = \tilde{q} + \Delta t g(\tilde{q})$$  这个分离是一个一阶精确算法，可以泰勒展开验证一下：$$q^{n+1} = (q^{n} + \Delta t f(q^{n})) + \Delta t g(q^{n} + \Delta t f(q^{n})) = q^n + \frac{dq}{dt}\Delta t + O(\Delta t^2)$$  正如你看到的这样，原来的微分方程被分离成两个单独的方程，一个只包含f函数，一个只包含g函数。需要注意的是 分离后的方程的解决顺序也有讲究，一般顺序解决，因为有的积分器需要另一个积分器的输入，顺序序列能保证不出错， 并行计算可能发生问题。2、N-S方程分解  有了对以上对$Navier-Stokes$方程的理论支撑，接下来我们就要如何用计算机来对该组偏微分方程进行离散化求解。为了程序的松耦合性以及使计算尽可能地高效、简单，在流体模型领域，我们将流体方程分成几个独立的步骤，然后按顺序先后推进。对于不可压缩的无粘度流体方程（即前面的欧拉方程$(2.19)$和$(2.20)$，我们将其离散化成对流项（advection）如公式$(3.1)$、体积力项（body force）如公式$(3.2)$、压力/不可压缩项如公式$(3.3)$：$$\frac{Dq}{Dt}=0 \tag {3.1}$$$$\frac{\partial \vec u}{\partial t}=\vec g \tag {3.2}$$$$ \begin{cases} \frac{\partial \vec u}{\partial t}+\frac{1}{\rho}∇p=0\\ ∇\cdot\vec u=0 \end{cases} \tag {3.3}$$  需要注意的是，在对流项公式$(3.1)$中我们用了一个通用量的符号$q$是因为我们不仅仅要对流体的速度进行对流，还需要对其他物理量进行对流。我们记对流项公式$(3.1)$的对流计算算法为$advect(\vec u, \Delta t, q)$，即对于给定的时间步长$\Delta t$和速度场$\vec u$，对物理量q进行对流。  对于体积力项$(3.2)$，我们采用简单的前向欧拉法即可：$\vec u \leftarrow \vec u + g\Delta t$。  对于压力/不可压缩项$(3.3)$，我们用一个称为$project(\Delta t, \vec u)$的算法，通过$project(\Delta t, \vec u)$计算出正确的压力以确保速度场$\vec u$的无散度性质。欧拉方案不会着重研究具体粒子间的作用力，因而不会正向去求解$\frac{1}{\rho}∇p$，它是利用流体不可压缩的特性，将速度场$\vec u$投影到散度为$0$的空间上，间接地解算了压力项。这种思想相当于，已知一个中间量$\vec u_{temp}$，对这个中间量的唯一一个操作（如正向求解压力$\frac{1}{\rho}∇p$）不可行，但是直到最终量$\vec u_{fianl}$符号的一个性质（散度为$0$），于是只要将$\vec u_{temp}$投影到符合散度为$0$的特性平面上，即可间接地还原正向求解压力的操作，得到最终的速度场$\vec u_{temp}$。  对流项$advect(\vec u, \Delta t, q)$的输入速度场$\vec u$要确保为无散度的状态，投影项$project(\Delta t, \vec u)$确保了流体体积保持不变，因而投影项输出的速度场必然是无散度的。所以我们只要确保投影项$project(\Delta t, \vec u)$输出的速度场$\vec u$作为对流项$advect(\vec u, \Delta t, q)$的输入即可，这时我们的分步求解流体方程的优势就体现出来了，其伪代码如下所示。算法1 Fluid Simulation($\vec u_n$, $\Delta t$):1: 初始化速度场$\vec u_n$,使得$\vec u_n$无散度 2: 对于每个时间步$n = 0,1,2,…$3:   决定一个合理的时间步长$\Delta t = t_{n+1}-t_n$ 4:   对流项计算$\vec u_A=advect(\vec u_n,\Delta t,\vec q)$ 5:   体积力项计算$\vec u_B=\vec u_A+\Delta t\vec g$6:   无散度投影$\vec u_{n+1}=project(\Delta t,\vec u_B)$3、时间步长  在流体模拟算法中，确定适当的时间步长是算法的第一步。因为计算流体模拟的最后结果是呈现在屏幕上的，所以$\Delta t$的选取与屏幕的刷新率有重要的关系。若选取的$\Delta t$有$t_n+\Delta t > t_{frame}$，那么必须做一个截断使$\Delta t=t_{frame}-t_n$。此外，流体模拟的三个步骤即对流项、体积力项、无散度投影项对时间步长$\Delta t$的要求不尽相同，要选择一个满足所有要求的最小时间步长能确保计算的收敛性。此外，一方面为了流体模拟的真实性，我们可能需要选取一个足够小的时间步长来复现流体的高质量细节。另一方面，有时高性能的需求又使得我们不能选取太小的时间步长去渲染一帧。假设一帧至少要进行三个时间步的模拟，那么$\Delta t$应该至少设成帧间隔时间的三分之一。4、网格结构  欧拉法的整个流程都是基于网格的，所以合理的网格结构是算法高效的关键点。$Harlow$和$Welch$提出了一种经典的$MAC$（marker and cell）网格结构，许多不可压缩流体模拟的算法都在这个网格结构上呈现出了良好的效率。$MAC$网格是一种交叉排列的网格，不同类型的物理量被存储于网格的不同位置。以二维的网格为例，如图3-1左图所示，流体粒子的压力数据存储于网格的中心点$P_{i,j}$，而速度则沿着笛卡尔坐标被分成了两部分。水平方向的$u$成分被存储在了网格单元竖直边的中心处，例如网格单元$(i,j)$和$(i+1,j)$之间的水平速度记为$u_{i+1/2,j}$。垂直方向的$v$成分则被存储在了网格单元水平面的中心上。这样的存储方案十分有利于估算流体流进/流出某个网格单元的量。图3-1 MAC网格,左图二维,右图三维  扩展到三维的情况，$MAC$网格同样是交错排列的结构网格，如图3-1右图所示。压力数值存储在立方体网格单元的中心，三个速度分量分别被记录在立方体网格单元的三个表面的中心点上。在数值计算时，这样的分配方式使得我们可以准确地采用中心差分法计算压力梯度和速度的散度，同时克服了中心差分法的一个普遍的缺点。一维的情况为例，在网格顶点位置$…,q_{i-1},q_i,q_{i+1}…$上估算量场$q$的导数，为了无偏（所谓无偏，就是不偏向左边或者右边）估计网格顶点$i$处的$\frac{\partial q}{\partial x}$，一种比较自然的方式就是采用一阶中心差分法：$$(\frac{\partial q}{\partial x})_i\approx \frac{q_{i+1}-q_{i-1}}{2\Delta x} \tag {3.4}$$  公式$(3.4)$是无偏的，且精确度为$O(\Delta x^2)$。而前向欧拉差分法偏向右边且精确度只有$O(\Delta x)$：$$(\frac{\partial q}{\partial x})_i\approx \frac{q_{i+1}-q_i}{\Delta x} \tag {3.5}$$  然而，公式$(3.4)$存在着一个非常严重的问题：网格点$i$的估算导数完全忽略了$q_i$的值。数学上，只有常数函数的一阶导数为零。但是公式$(3.4)$遇到了锯齿函数如$q_i=(-1)^i$时，它错误地将该类函数的导数估算为$0$，这种问题被称为零空间问题（null-space problem）。  交叉错排的$MAC$网格完美地克服了中心差分法的零空间问题，同时也保持了它的无偏二阶精度。在$MAC$网格上运用中心差分法，网格点$i$处的估算导数公式如下所示：$$(\frac{\partial q}{\partial x})_i\approx\frac{q_{i+1/2}-q_{i-1/2}}{\Delta x} \tag {3.6}$$  $MAC$网格确实给流体的压力计算和不可压缩性的处理带来了很大的便利，但与此同时也带来了一些其他方面的麻烦。如果我们要估算某个地方的速度向量，即便采样点恰好在网格点上我们也要做一些插值才能获取相应的速度向量。在网格点处，我们通常采用平均法，以二维为例：$$\vec u_{i,j}=(\frac{u_{i-1/2,j}+u_{i+1/2,j}}{2},\frac{v_{i,j-1/2}+v_{i,j+1/2}}{2}),\\ \vec u_{i+1/2,j}=(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}),\\ \vec u_{i,j+1/2}=(\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}).\tag {3.7}$$  最后，在实现中下标索引一般没有浮点数之说，前面直接采用$i+1/2$的记法是为了便于叙述。一般约定如下：$$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}. \tag{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)$数组中。四、对流算法  求解如下所示的对流方程是流体模拟的关键一步：$$\frac{Dq}{Dt}=0 \tag {4.1}$$  我们把这个对流数值计算的算法记为：$$q^{n+1}=advect(\vec u,\Delta t,q^n) \tag {4.2}$$  公式$(4.2)$中的各个符号含义：  $\vec u$：在$MAC$网格上的离散化的速度场；  $\Delta t$：时间步长；  $q^n$：当前的物理量场$q$（如流体密度、速度、燃烧物浓度等）；  $q^{n+1}$：经过对流后得到的新的量场。  在这里要特别注意，输入对流算法的速度场$\vec u$必须是无散度的，否则模拟结果会出现一些奇怪的失真现象。1、半拉格朗日对流算法（Semi-Lagrangian Advection）  一维情况下，对流方程$(4.1)$写成偏微分的形式如下：$$\frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=0 \tag {4.3}$$  分别采用前向欧拉差分法计算对时间的偏导和中心差分法计算对空间的偏导，我们有：$$\frac{q^{n+1}_{i}-q^n_i}{\Delta t}+u^n_i\frac{q^n_{i+1}-q^n_{i-1}}{2\Delta x}=0 \tag {4.4}$$  转成以$q^{n+1}_i$为计算目标的显式公式，得：$$q^{n+1}_i=q^n_i-\Delta t u^n_i\frac{q^n_{i+1}-q^n_{i-1}}{2\Delta x} \tag {4.5}$$  公式$(4.5)$看起来没什么问题，但是却存在非常严重的漏洞。首先，前向欧拉法被证明是无条件不稳定的空间离散方法：无论取多么小Δ𝑡，随着时间步的推进，累积误差终将发散。即使使用更稳定的时间积分方法来取代前向欧拉方法，解决了时间上的PDE（Partial Differential Equation，偏微分方程）计算，空间上的PDE计算还是会带来重大的麻烦。标准中心差分方法不可避免地会出现的零空间问题，具有高频震荡性质的速度场对空间的导数被错误地计算为$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$个时间步时$\vec x_G$的$q^{n+1}_G$即为第$n$个时间步时$\vec x_P$的$q^{n}_P$。如下图4-1为半拉格朗日对流法的示意图。图4-1 半拉格朗日对流法  半拉格朗日对流法的第一步就是要找出$\vec x_P$，为此我们根据$\vec x_G$做反向的追踪。粒子位置对时间的导数就是速度场：$$\frac{d\vec x}{dt}=\vec u(\vec x) \tag {4.6}$$  经过一个时间步长$\Delta t$之后，粒子由$\vec x_P$移动到$\vec x_G$。为了得到$\vec x_P$，最简单的方法就是采用前向欧拉法进行倒推：$$\vec x_P=\vec x_G-\Delta t\vec u(\vec x_G) \tag {4.7}$$  然而前向欧拉法只有一阶的精度，若在不改变$\Delta t$的情况下提高精度，我们可以采用高阶的龙格库塔法（Runge-Kutta method）。采用二阶的龙格库塔法如下所示：$$\vec x_{mid}=\vec x_G-\frac12\Delta t\vec u(\vec x_G),\\ \vec x_P=\vec x_G-\Delta t\vec u(\vec x_{mid}). \tag {4.7}$$  倒推得到$\Delta t$之前的网格位置$\vec x_P$一般不会恰好在网格顶点上，为此我们需要做些插值。三维模拟通常采用三线性插值，而二维的则采用双线性插值。$$q^{n+1}_G=interpolate(q_n,\vec x_P) \tag {4.8}$$2、边界情况  若我们倒推得到的$\vec x_P$仍然在流体的内部，那么做插值是完全没问题的。但若$\vec x_P$在流体的边界之外呢？这种情况的出现的原因通常有两个：一个是$\vec x_P$确确实实在流体的外部且即将流入流体内部，另一个是由前向欧拉法或龙格库塔法的数值计算方法带来的误差导致。  在一种情况下，我们应该知道当流体流入时其携带的物理量，此时我们将这个外部流入的物理量作为返回值即可。例如，第$n$个时间步时的外部流体以速度$\vec U$和温度$T$在第$n+1$个时间步时注入流体内部$\vec x_G$的位置，那么$\vec T^{n+1}_G$的值就为$T$。  在第二种由误差导致的情况下，一个适当的策略就是根据边界上的最近点外推出所求得物理量。在模拟某些流体时，外推变得很简单。例如，在模拟烟雾时我们简单地假设烟雾流体外部即空气的速度风场为某个常数$\vec U$（可能为$0$），这样边界上的速度场都取$\vec U$。但还有一些必须根据流体内部的已知量外推出未知量，这时情况就变得比较复杂了。具体如何外推将在后面介绍，目前我们只需要知道大概的步骤：首先寻找边界上的最近点，然后在最近点的领域内插值获取相应的物理量场。3、时间步长大小  对任何一种数值计算方法的主要的考虑点就是它是否稳定。幸运的是，半拉格朗日对流法已经被证明是一种无条件稳定的算法：无论$\Delta t$取多大，它永远不会出现数值爆炸的现象。因为每一个新值$q$的确定，都是通过对旧值得插值，无论是线性插值、双线性插值还是三线性插值，$q$的大小都是处于插值点之间，不会得到比原来插值点更大或者更小的值，因而$q$是有上下界的。这使得我们可以尽情地根据所需的模拟质量和模拟效率去调整时间步长。  但是在实践中，时间步长的大小也不能选得太过极端，否则会产生一些奇观的现象。Foster和Fekiw提出了一个对$\Delta t$的限制：流体粒子在$\Delta t$内的倒推轨迹最多经过某个常数个网格单元为宜，例如5个：$$\Delta t \leq \frac{5\Delta x}{u_{max}} \tag {4.9}$$  公式$(4.9)$中，$u_{max}$是速度场的最大值，我们可以简单地取 存储在网格中的最大速度值。一个更鲁棒的方法考虑了体积力（如重力、浮力等）对最大速度的影响：$$u_{max}=max(|u^n|)+\Delta t|g| \tag {4.10}$$  将不等式$(4.9)$的最大值带入公式$(4.10)$，我们有：$$u_{max}=max(|u^n|)+\frac{5\Delta x}{u_{max}}|g| \tag {4.11}$$  取一个简单的速度上界（简化了公式$(4.11)$），$u_{max}$：$$u_{max}=max(|u^n|)+\sqrt{5\Delta xg} \tag {4.12}$$  这样确保了$u_{max}$始终为正，且避免公式$(4.9)$的除$0$错误。  关于时间步长的讨论离不开$CFL$（以Courant、Friedrichs、Lewy三人的名字命名）条件。$CFL$条件是一个简单而直观的判断计算是否收敛的必要条件。它的直观物理解释就是时间推进求解的速度必须大于物理扰动传播的速度，只有这样才能将物理上所有的扰动俘获到。满足$CFL$条件意味着当$\Delta x$和$\Delta t$趋于取极限$0$时，数值计算所求的解就会收敛到原微分方程的解。  对于半拉格朗日对流法，其满足$CFL$条件当且仅当在极限情况下，追踪得到的粒子轨迹足够逼近真实的轨迹。足够逼近的意思是经过正确的网格插值能够得到正确的依赖域（即差分格式的依赖域包含了原微分方程的依赖域），追踪的轨迹就会收敛到正确真实的轨迹。  因而，对于采用标准的显式有限差分法的对流方程求解，为了保证收敛，我们要求$q^{n+1}$的新值是由以当前网格点为中心、以$C\Delta x$（$C$是一个小的整数常量）为半径的邻域范围内插值得到：$$\Delta t \leq C\frac{\Delta x}{|\vec u|} \tag {4.13}$$  公式$(4.13)$中的$C$被称为$CFL$数，因而不等式$(4.9)$可以看成是公式$(4.13)$取$CFL$数为$5$得到。4、数值耗散  对流算法在对流获取新的物理量场$q^{n+1}_i$时会进行一些插值操作，插值不可避免地会平滑物理量场，这带来了一些数值耗散。一次两次的数值耗散不会由太大的影响，但是在流体模拟中我们会在每个时间步都进行对流运算，反反复复的平滑操作将数值耗散不断扩大，损失大量的流体细节。  以一维的对流项计算为例，流体速度为常量$u>0$：$$\frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=0 \tag {4.14}$$  假设$\Delta t < \frac{\Delta x}{u}$，即单个时间步长内粒子追踪轨迹长度小于单个网格单元的大小。我们的目标点是$x_i$，则倒推得到的粒子位置就落在了$[x_{i-1},x_i]$上的$x_i-\Delta tu$，然后进行线性插值得到$q^{n+1}_i$：$$q^{n+1}=\frac{\Delta tu}{\Delta x}q^n_{i-1}+(1-\frac{\Delta tu}{\Delta x})q^n_i \tag {4.15}$$  将公式$(4.15)$整理一下，有：$$q^{n+1}_i=q^n_i-\Delta tu\frac{q^n_i-q^n_{i-1}}{\Delta x} \tag {4.16}$$  公式$(4.16)$实际上正好就是采用时间上的前向欧拉差分法和空间上的单向有限差分法的欧拉方案，把$q^n_i$看成是$q^n$关于$x_i$的函数，对$q^n_{i-1}$进行泰勒级数展开：$$q^n_{i-1}=q^n_i-(\frac{\partial q}{\partial x})^n_i\Delta x+(\frac{\partial^2q}{\partial x^2})^n_i\frac{\Delta x^2}{2}+O(\Delta x^3) \tag {4.17}$$  将公式$(4.17)$代入公式$(4.16)$，并做一些变量消去，可得：$$ q^{n+1}_i=q^n_i-\Delta tu(\frac{\partial q}{\partial x})^n_i+\Delta tu\Delta x(\frac{\partial^2q}{\partial x^2})^n_i+O(\Delta x^2) \tag {4.18} $$  在二阶截断误差的情况下，结合公式$(4.18)$和公式$(4.14)$，有：$$\frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=u\Delta x(\frac{\partial^2q}{\partial x^2}) \tag {4.19}$$  右边就是对流方程计算时引入的额外类似粘度乘上系数$u\Delta x$的项。这也就是说，当我们采用简单的半拉格朗日法去求解无粘度的对流方程时，模拟的结果却看起来我们像时在模拟有粘度的流体。 这就是数值耗散！当然，当$\Delta x\to 0$时，这个数值耗散系数也会趋于$0$，所以取时间步无穷小时能够得到正确的模拟结果，但这需要耗费巨额的计算资源开销。我们通常模拟的流体大多数都是无粘度的，所以如何减少这个数值耗散是个至关重要的难题。  一个简单有效的修复数值耗散的方法就是采用更加锐利的插值方法，从而尽可能地减少由插值带来的数值耗散。在一维的情况时，我们采用三次插值（cubic interpolant）如下公式$(4.21)$，而不是简单的一次线性插值$(4.20)$：$$q\approx(1-s)x_i+sx_{i+1} \tag {4.20}$$$$q\approx[-\frac13s+\frac12s^2-\frac16s^3]q_{i-1}+[1-s^2+\frac12(s^3-s)]q_i\\ +[s+\frac12(s^2-s^3)]q_{i+1}+[\frac16(s^3-s)]q_{i+2} \tag {4.21}$$  扩展到二维或者三维就是双三次插值（bicubic interpolation）或三三次插值（tricubic interpolation）。以二维情况为例，我们可以先沿着$x$轴做第一遍的三次插值如公式$(4.22)$，然后再沿着$y$轴做第二遍插值如公式$(4.23)$：$$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}=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}. \tag {4.22}$$$$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}$$  当然也可以先沿着$y$轴，然后再沿着$x$轴做插值操作。五、求解流体不可压缩的泊松方程在进行半拉格朗日对流之后，还需要对流体速度进行修正，使流体的速度场散度为零。主要是关于流体压强的泊松方程求解，通常采用共轭梯度法求解大规模稀疏矩阵的线性方程组，辅以不完全的Cholesky因子化。涉及到的数学内容比较多。离散的压力梯度离散的速度场散度求解关于压强的泊松方程不可压缩投影参考资料1、流体的不可压缩投影  在流体模拟中，确保流体的不可压缩性是非常重要、关键的一步，也是最耗时的一步，这是整个流体模拟过程中的核心，涉及到模拟效果的真实性、模拟过程的效率快慢两方面的内容，因而相关的数学内容比较多。在流体模拟中，确保流体的不可压缩性的算法我们通常用$project(\Delta t, \vec u)$来表示，它输入时间步长$\Delta t$和流体的速度场$\vec u$，然后在流体内部的压力的作用下更新速度场：$$\vec u^{n+1}=\vec u-\Delta t\frac1\rho \nabla p \tag {5.1}$$  上述公式中，$\vec u^{n+1}$就是新的速度场，$\rho$是流体密度，$p$为压力。通过公式$(5.1)$，我们就可以得到满足流体不可压缩性的流体，此时得到的新速度场$\vec u^{n+1}$的散度为0：$$\nabla \cdot \vec u^{n+1}=0 \tag {5.2}$$  且在固体边界处，其速度场在法线方向上的投影等于固定墙的速度场在法线方向上的投影：$$\vec u^{n+1}\cdot \vec n=\vec u_{solid}\cdot \vec n \tag {5.3}$$  在前一步展开如何求解上述的几个公式之前，我们先把流体模拟区域划分成一个一个格子，有流体的格子标记为F，固体边界区域的格子标记为$S$，空白的区域则标记为$E$。下图1展示了一个二维的例子。图1 模拟区域标记   模拟的网格依旧是MAC网格，即交错的网格结构，如下图2所示。格子中心存储流体压力值，边缘部分存储流体的速度场向量，依次交错存储。图2 MAC网格2、离散的压力梯度  MAC网格解决了普通的网格的零域（null space）问题，在使用中心差分法计算梯度时鲁棒性更强。目前暂时假设流体的压力值已经知道且存储在MAC网格上，那么二维和三维的压力梯度计算采用中心差分法如下所示：$$\nabla p=(\frac{p_{i+1,j}-p_{i,j}}{\Delta x},\frac{p_{i,j+1}-p_{i,j}}{\Delta x})\\ \nabla p=(\frac{p_{i+1,j,k}-p_{i,j,k}}{\Delta x}, \frac{p_{i,j+1,k}-p_{i,j,k}}{\Delta x}, \frac{p_{i,j,k+1}-p_{i,j,k}}{\Delta x}) \tag {5.4}$$  然后将公式$(5.4)$中的离散压力梯度公式代入公式$(5.1)$，可以得到实际上的速度场更新公式：$$\begin{cases} u_{i+1/2,j}^{n+1}=u_{i+1/2,j}-\Delta t\frac1\rho\frac{p_{i+1,j}-p_{i,j}}{\Delta x}\\ v_{i,j+1/2}^{n+1}=v_{i,j+1/2}-\Delta t\frac1\rho\frac{p_{i,j+1}-p_{i,j}}{\Delta x} \end{cases}\\ \begin{cases} u_{i+1/2,j,k}^{n+1}=u_{i+1/2,j,k}-\Delta t\frac1\rho\frac{p_{i+1,j,k}-p_{i,j,k}}{\Delta x}\\ v_{i,j+1/2,k}^{n+1}=v_{i,j+1/2,k}-\Delta t\frac1\rho\frac{p_{i,j+1/2,k}-p_{i,j,k}}{\Delta x}\\ w_{i,j,k+1/2}^{n+1}=w_{i,j,k+1/2}-\Delta t\frac1\rho\frac{p_{i,j,k+1/2}-p_{i,j,k}}{\Delta x} \end{cases} \tag {5.5}$$  公式$(5.5)$给出了二维和三维的流体不可压缩速度场更新公式。注意，固体墙边界并没有压力一说，而自由面边界压力我们设置为零，这样上述的公式仅适用于那些不邻近固体边界的速度场分量，一般情况下至少需要邻近一个流体区域。  现在我们把目标放到边界条件的处理上。首先是自由面，这部分通常就是空气，我们直接设置其压力值为零，这类边界处理我们称之为第一类边界条件——狄利克雷边界条件（Dirichlet boundary condition）$^{[1]}$。所谓狄利克雷边界条件，就是我们直接设定边界处的值，这是在数值方法中最简单的一种边界处理方法。  然后就是固体墙边界的处理。在固体墙边界部分，我们需要将流体的速度场在法线方向上的投影等于固体边界在法相方向上的速度场向量值，一种方法就是直接设置，例如在$(i+1/2,j,k)$处的速度场分量$u$，假设其左边$(i,j,k)$为固体边界区域，右边$(i+1,j,k)$为流体区域，那么直接赋值$u$：$$u_{i+1/2,j,k}^{n+1}=u_{i+1/2,j,k}^{solid} \tag {5.6}$$  直接赋值需要一个额外的分支进行。另外一种方案就是给固定墙区域设定一个ghost压强值，只要通过公式$(5.5)$更新后固体边界处的速度场为公式$(5.6)$即可，此时不需要额外的分支。通过一些变换，可得固体墙区域的ghost压强计算公式：$$p_{i,j,k}^{ghost}=p_{i+1,j,k}-\frac{\rho\Delta x}{\Delta t}(u_{i+1/2,j,k}-u_{i+1/2,j,k}^{solid}) \tag {5.7}$$  这种方案在进行不可压缩投影时不需要考虑邻域是否是固体边界，统一用公式$(5.5)$计算，在固体边界处用ghost压强值计算，依旧以上面为例：  实际上，将公式$(5.7)$代入公式$(5.8)$中，就可以得到公式$(5. 6)$，通过这种方式不再将固体边界的处理作为一个单独考虑的情况。将公式$(5.8)$做一些移项，重写为：$$\frac{\Delta t}{\rho}\frac{p_{i+1,j,k}-p_{i,j,k}^{ghost}}{\Delta x}= u_{i+1/2,j,k}-u_{i+1/2,j,k}^{solid} \tag {5.9}$$  上述公式中的左边压强部分可以看成是关于压强梯度的有限差分，即关于下面公式的一个近似：$$\frac{\Delta t}{\rho}\frac{\partial p}{\partial x}=u-u^{solid} \tag {5.10}$$  公式$\frac{\partial p}{\partial x}$仅仅是关于一个方向，同样地，我们可以取其他方向的固体边界条件，得到更一般的公式：$$\frac{\Delta t}{\rho}\nabla p\cdot \vec n=(\vec u-\vec u^{solid})\cdot \vec n \tag {5.11}$$  公式$(5.11)$给出了固体边界处的压力梯度值，与狄利克雷条件处理方式不同，这里给出的是梯度的取值，这种边界处理方式被称为第二类边界条件——冯诺依曼边界条件（Neumannn boundary condition）$^{[1]}$。最后，假设流体的压力已知（实际上未知），根据压力梯度更新流体的速度场，算法伪代码如下：scale = dt/(density*dx); loop over i,j,k: #update u component. if label(i-1,j,k)==FLUID or label(i,j,k)==FLUID: if label(i-1,j,k)==SOLID or label(i,j,k)==SOLID: u(i,j,k) = usolid(i,j,k); else u(i,j,k) -= scale * (p(i,j,k) - p(i-1,j,k)); else mark u(i,j,k) as unknown; #update v component. if label(i,j-1,k)==FLUID or label(i,j,k)==FLUID: if label(i,j-1,k)==SOLID or label(i,j,k)==SOLID: v(i,j,k) = vsolid(i,j,k); else v(i,j,k) -= scale * (p(i,j,k) - p(i,j-1,k)); else mark v(i,j,k) as unknown; #update w component. if label(i,j,k-1)==FLUID or label(i,j,k)==FLUID: if label(i,j,k-1)==SOLID or label(i,j,k)==SOLID: w(i,j,k) = wsolid(i,j,k); else w(i,j,k) -= scale * (p(i,j,k) - p(i,j,k-1)); else mark w(i,j,k) as unknown;3、离散的速度场散度  流体不可压缩的最终表现就是速度场的散度为零，即$\nabla \cdot \vec u=0$。经过不可压缩投影之后新的速度场$\vec u^{n+1}$应该满足散度为零的条件，我们通过有限差分法估算新速度场$\vec u^{n+1}$的散度。首先是散度的数学定义，二维和三维的散度定义如下所示：$$\nabla \cdot \vec u=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\\ \nabla \cdot \vec u=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+ \frac{\partial w}{\partial z} \tag {5.12}$$  给定一个二维欧拉网格$(i,j)$，采用中心差分法计算离散的散度公式如下：$$(\nabla\cdot \vec u)_{i,j}\approx \frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta x}+\frac{v_{i,j+1/2}-v_{i,j-1/2}}{\Delta x} \tag {5.13}$$  同理，三维$(i,j,k)$的离散散度计算公式为：$$(\nabla\cdot \vec u)_{i,j,k}\approx \frac{u_{i+1/2,j,k}-u_{i-1/2,j,k}}{\Delta x} +\frac{v_{i,j+1/2,k}-v_{i,j-1/2,k}}{\Delta x} +\frac{w_{i,j,k+1/2}-w_{i,j,k-1/2}}{\Delta x} \tag {5.14}$$  事实上，我们更关心散度的相反数，且仅对于流体区域计算其散度值，其余的固体边界、自由面部分不需要计算散度值。下面是一个计算流体散度的伪代码。scale = 1/dx; loop over i,j,k where label (i,j,k)==FLUID: rhs(i,j,k)=-scale * (u(i+1,j,k)-u(i,j,k) +v(i,j+1,k)-v(i,j,k) +w(i,j,k+1)-w(i,j,k));  给定一个流体区域及其表面$\partial cell$，流体散度的意义就是流体流进流出的总比例：$$\int\int_{\partial cell}\vec u \cdot \vec n \tag {5.15}$$  若流体区域为一个欧拉网格的立方体格子，那么上述的积分就是对立方体六个面上的速度场进行积分。这种方法与我们前面讨论的有限差分法不同，这是有限体积法（Finite volume method）。这里只是顺带一提。4、求解关于压强的泊松方程  在前面我们讨论了如何在已知流体内部压强的情况下进行不可压缩投影，更新速度场，并估算速度场的散度值。实际上，流体内部压强并不知道，这需要我们计算得到。接下来的内容就是关于如何求解流体的压力数值，这部分涉及的数学内容比较多，是比较难的一部分。泊松方程  我们知道新的速度场$\vec u^{n+1}$需要满足散度为零，从而保证流体的不可压缩条件。又已知新的速度场$\vec u^{n+1}$计算公式为上面的公式$(5.5)$，离散的散度计算公式为公式$(5.13)$和公式$(5.14)$，我们将公式$(5.5)$得到的速度场代入公式$(5. 13)$、$(5.14)$中，就可以得到以流体压力为未知量的线性方程。  我们先来看二维的情况，关于$\vec u^{n+1}$的散度为零表示成如下：$$\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 \tag {5.16}$$  将公式$(5.5)$中二维的速度更新公式代入上面的公式$(5.16)$中，可得：$$\frac{1}{\Delta x}[(u_{i+1/2,j}-\Delta t\frac1\rho\frac{p_{i+1,j}-p_{i,j}}{\Delta x}) -(u_{i-1/2,j}-\Delta t\frac1\rho\frac{p_{i,j}-p_{i-1,j}}{\Delta x})\\ +(v_{i,j+1/2}-\Delta t\frac1\rho\frac{p_{i,j+1}-p_{i,j}}{\Delta x}) -(v_{i,j-1/2}-\Delta t\frac1\rho\frac{p_{i,j}-p_{i,j-1}}{\Delta x})]=0\\ \to \\ \frac{\Delta t}{\rho}(\frac{4p_{i,j}-p_{i+1,j}-p_{i,j+1}-p_{i-1,j}-p_{i,j-1}}{\Delta x^2}) =-(\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta x}+\frac{v_{i,j+1/2}-v_{i,j-1/2}}{\Delta x}) \tag {5.17}$$  三维的情况同理，只不过多了一维需要处理：$$\frac{u_{i+1/2,j,k}^{n+1}-u_{i-1/2,j,k}^{n+1}}{\Delta x} + \frac{v_{i,j+1/2,k}^{n+1}-v_{i,j-1/2,k}^{n+1}}{\Delta x} + \frac{w_{i,j,k+1/2}^{n+1}-w_{i,j,k-1/2}^{n+1}}{\Delta x} = 0 \tag {5.18}$$$$\frac{1}{\Delta x}[(u_{i+1/2,j,k} -\Delta t\frac1\rho\frac{p_{i+1,j,k}-p_{i,j,k}}{\Delta x}) -(u_{i-1/2,j,k}-\Delta t\frac1\rho\frac{p_{i,j,k}-p_{i-1,j,k}}{\Delta x})\\ +(v_{i,j+1/2,k}-\Delta t\frac1\rho\frac{p_{i,j+1,k}-p_{i,j,k}}{\Delta x}) -(v_{i,j-1/2,k}-\Delta t\frac1\rho\frac{p_{i,j,k}-p_{i,j-1,k}}{\Delta x})]\\ +(w_{i,j,k+1/2}-\Delta t\frac1\rho\frac{p_{i,j,k+1}-p_{i,j,k}}{\Delta x}) -(w_{i,j,k-1/2}-\Delta t\frac1\rho\frac{p_{i,j,k}-p_{i,j,k-1}}{\Delta x})]=0\\ \to \\ \frac{\Delta t}{\rho}(\frac{6p_{i,j,k}-p_{i+1,j,k}-p_{i,j+1,k}-p_{i,j,k+1}-p_{i-1,j,k}-p_{i,j-1,k}-p_{i,j,k-1}}{\Delta x^2}) \\ = -( \frac{u_{i+1/2,j,k}-u_{i-1/2,j,k}}{\Delta x} + \frac{v_{i,j+1/2,k}-v_{i,j-1/2,k}}{\Delta x} + \frac{w_{i,j,k+1/2}-w_{i,j,k-1/2}}{\Delta x} ) \tag {5.19}$$  仔细观察上面的公式$(5.17)$和公式$(5.19)$，可以发现他们都是泊松问题$-\Delta t/\rho \nabla \cdot \nabla p=-\Delta \cdot \vec u$的一个离散计算公式，左边是关于压力的拉普拉斯算子的数值近似，右边是关于速度场散度的数值近似，我们最终要求解的就是这么一个关于压力的方程。方程右边是散度的相反数，这就是我们前面计算并保存的是散度的相反数的原因。当邻居格子是自由面边界时，我们直接将其相应的压力值取零；当邻居格子是固体墙边界时，我们根据冯诺依曼边界条件计算得到的压强替换相应的压力项。写成矩阵形式  上面讨论的仅仅是一个格子的压力项求解方程，事实上我们要求解的是整个流体区域的压力项值。为此，为了便于阐述和求解，我们将上面的线性方程写成矩阵乘向量的形式。我们将方程中所有压力项的系数提出来，写成一个系数矩阵$A$，所有流体区域的压力项写成一个未知变量的向量$p$，相应的所有流体区域的负散度作为线性方程组右边的向量$b$，那么我们要求解的就是下面的一个大规模非齐次线性方程组：$$Ap=b \tag {5.20}$$  矩阵$A$每一行存储的是一个流体格子求解压力项的系数。仔细观察公式$(5.17)$和$(5. 19)$，实际上每一行的系数中仅仅只有几个不为0，不为0的系数是其周围邻居和自身的。在三维情况下，给定$(i,j,k)$这个流体格子，那么其对应的矩阵$A$中的那一行中，仅仅是$p_{i,j,k}$、$p_{i\pm 1,j,k}$、$p_{i,j\pm 1,k}$和$p_{i,j,k\pm 1}$对应的系数不为0，即又7个不为的系数，剩下的系数全部为0。二维就是一行有5个系数不为0。可以看到矩阵$A$是一个大规模的稀疏矩阵，大部分的矩阵元素都为0，所以我们没有必要直接存储矩阵$A$，因为这样将会极大地耗费内存，甚至出现内存不足的情况，因为矩阵A实在是太大了。   仔细观察前面公式$(5. 17)$和$(5.19)$的线性方程，对于给定的流体格子$(i,j,k)$，其邻域压力项的系数为$-\Delta t/(\rho \Delta x^2)$，自身压力项的系数为邻居格子中流体格子或者空气（自由面）格子的数量$n_{i,j,k}$在乘上$\Delta t/(\rho\Delta x^2)$，即$n_{i,j,k}\Delta t/(\rho \Delta x^2)$，可以看到在三维时$n_{i,j,k}$至多为6、在二维时至多为$4$。  同时，系数矩阵$A$也是一个对阵矩阵。举个例子，$A_{(i,j,k)(i+1,j,k)}$是方程组中关于$p_{i+1,j,k}$的系数，那么它必然等于$A_{(i+1,j,k)(i,j,k)}$，这是因为$(i+1,j,k)$是$(i,j,k)$的邻居，那么$(i,j,k)$必然也是$(i+1,j,k)$的邻居。系数矩阵$A$既然是对阵矩阵，那么我们存储系数矩阵$A$时又可以减少一半的存储量了，只需存储上三角形、或者下三角矩阵。  经过上述的讨论，对于大规模的稀疏矩阵$A$，以二维为例，我们采用的存储方案为：每一个流体格子存储关于自身压力项的系数$A_{(i,j)(i,j)}$，这个就是矩阵$A$中的对角线元素，而关于邻居流体格子压力项的系数我们仅存储正向方向上的（因为对称）$A_{(i,j),(i+1,j)}$和$A_{(i,j)(i,j+1)}$。我们记矩阵$A$的对角线元素为Adiag(i,j)，$x$轴方向的矩阵元素Ax(i,j)，$y$轴方向的矩阵元素Ay(i,j)。三维同理，分别存储到Adiag(i,j,k)、Ax(i,j,k)、Ay(i,j,k)、Az(i,j,k)，一个三维的系数矩阵$A$的计算并存储的伪代码如下所示：scale = dt / (density*dx*dx); loop over i,j,k: if label(i,j,k)==FLUID: # handle negative x neighbor if label(i-1,j,k)==FLUID: Adiag(i,j,k) += scale; # handle positive x neighbor if label(i+1,j,k)==FLUID: Adiag(i,j,k) += scale; Ax(i,j,k) = -scale; else if label(i+1,j,k)==EMPTY: Adiag(i,j,k) += scale; # handle negative y neighbor if label(i,j-1,k)==FLUID: Adiag(i,j,k) += scale; # handle positive y neighbor if label(i,j+1,k)==FLUID: Adiag(i,j,k) += scale; Ay(i,j,k) = -scale; else if label(i,j+1,k)==EMPTY: Adiag(i,j,k) += scale; # handle negative z neighbor if label(i,j,k-1)==FLUID: Adiag(i,j,k) += scale; # handle positive z neighbor if label(i,j,k+1)==FLUID: Adiag(i,j,k) += scale; Az(i,j,k) = -scale; else if label(i,j,k+1)==EMPTY: Adiag(i,j,k) += scale;共轭梯度算法   现在我们把目标放到求解大规模的非齐次线性方程组上，即上面的公式$(5.20)$。系数矩阵$A$除了是一个大规模的稀疏对称矩阵之外，还是一个正定（Positive definite）矩阵，即对于任意的非零向量$q$，$q^TAq>0$均成立，且系数矩阵$A$的所有特征值均大于0。当然在实际情况下，系数矩阵$A$可能不是一个严格的正定矩阵，而是一个半正定矩阵，即对于任意非零向量$q$，有$q^TAq\geq0$。这种情况出现在流体模拟区域中，存在一些流体区域完全被固体墙边界包围，此时的系数矩阵$A$不是一个严格的正定矩阵，甚至此时的矩阵$A$是一个奇异矩阵，它不存在逆矩阵，从而导致方程$(5. 20)$不一定有解。因而在一些特殊情况下，我们求解的是兼容条件下的解，即只要使得流体-固体边界处流体的流进和流出保持平衡（流进量等于流出量）即可。    求解一个线性方程组最常用、最暴力的方法就是高斯消元法，逐行进行消元，最后使得矩阵变成一个上三角矩阵，从最后一行往回代即可求解出方程的解（只要解存在）。但是高斯消元法实在是太慢了，小规模的矩阵来说还好，大规模的矩阵方程求解一般不会这么直接暴力解。这里我们采用的求解方程$(5.20)$的方法是共轭梯度算法（Conjugate gradient method，简称为CG）$^{[2]}$。相比于普通的高斯消元法，共轭梯度算法采用了一个迭代的流程，迭代到给定的收敛程度即可，整个算法流程只涉及到矩阵-向量乘法、向量加减、向量-标量乘法以及非常少的点乘操作，这些都可以快速并行地计算。    共轭梯度算法就是针对对称正定矩阵的线性方程组的求解，给定一个线性方程组如下$(5.21)$所示，其中系数矩阵$A$是一个元素均已知的实对称正定矩阵，右边的向量$\vec b$也已知，我们要求解的就是未知向量$\vec x$。这个就是共轭梯度所要解决的线性方程组。$$A x= b \tag {5.21}$$  给定两个非零向量$u$和$v$，我们说它们是共轭的（相对于$A$），若：$$u^TAv = 0 \tag {5.22}$$  而矩阵$A$是一个对称的正定矩阵（故有$A^T=A$），那么上式的左边可以写成多种向量内积的形式，它们完全是等价的。这样两个非零向量是共轭的，当且仅当这两个向量相对于下面的内积是正交的。共轭是一个对称的相互关系，$u$共轭于$v$，那么$v$也共轭于$u$。$$u^TAv=<u,v>_A=<Au,v>=<u,A^Tv>=<u,Av> \tag{5.23}$$  令$P=\{p_1,…,p_n\}$是一个这样的向量集合，集合中的所有向量都相对于$A$互相共轭，那么$P$构成了$R^n$空间的一个基，方程$(5.21)$的解$x_*$可以用这组基的线性组合来表示：$$x_*=\Sigma_{i=1}^n\alpha_i p_i \tag {5.24}$$  将上述的$x_*$表达式代入方程$(5.21)$中，有：$$Ax_*=\Sigma_{i=1}^n\alpha_iAp_i \tag {5.25}$$  再左乘上$P$中的一个向量$p_k^T$：$$p_k^TAx_*=\Sigma_{i=1}^n\alpha_ip_k^TAp_i \tag {5.26}$$  然后将$Ax_*=b$、$u^TAv=(u,v)_A$代入公式$(5.26)$中：$$p_k^Tb=\Sigma_{i=1}^n\alpha_i {<p_k,p_i>}_A \tag{5.27}$$  最后根据$u^Tv=(u,v)$以及$\forall i\neq k: (p_k,p_i)_A=0$，可得：$$ {<p_k,p_k>}=\alpha_k {<p_k,p_k>}_A \\ \to \\ \alpha_k=\frac{<p_k,b>}{<p_k,p_k>_A} \tag{5.28}$$  由此，我们可以得到一个求解方程$Ax=b$的方法，首先找到$n$个互相共轭的向量，然后计算其组合系数$\alpha_k$，最后未知向量就由这$n$个互相共轭的向量线性组合而成（即公式$(5.24)$）。但显然直接寻找$n$个互相共轭的向量在$n$比较大的时候不现实，因为这将耗费大量的时间。  共轭梯度采用的是迭代法，我们给出了一个初始预测解$x_0$，致力于寻找最终解$x_d$，通过迭代的方法使得初始的预测解$x_0$逐步向$x_d$逼近，每一次的迭代中我们挑选出一个共轭向量$p$。迭代过程中的目标就是最小化预测向量与最终解$x_d$之间的差距，此时我们的问题就转换为如下的二次方程的最优化问题：$$f(x)=x^T(\frac12Ax-b)=\frac12x^TAx-x^Tb, x\in R^n \tag {5.29}$$  公式$(5.29)$是一个关于$n$维向量的二次方程，是一个$n$元函数，其一阶导数和二阶导数分别如下所示：$$\nabla f(x)=Ax-b \\ \nabla^2 f(x)=A \tag {5.30}$$  可以看到当二次方程的一阶导数$Ax-b=0$时，此时的$x$就是我们要求的线性方程组的解，此时二次方程取得最小值。共轭梯度法在初次迭代时直接取第一个共轭向量$p_0=b-Ax_0$，即负梯度方向，这一步与梯度下降法类似。但是后续迭代取的向量就不再是负梯度向量了，而是取一个与梯度向量共轭的向量，这就是名称共轭梯度法的由来。实际上，每一次迭代时我们取的负梯度向量就是当前预测解与最终解之间的残差（Residual），第$k$次迭代时的残差为：$$r_k=b-Ax_k \tag {5.31}$$  为了确保在每一次迭代中取得的向量共轭于梯度以及前面迭代得到的向量，我们结合残差向量以及前面迭代的向量计算当前迭代挑选的向量：$$p_k=r_k-\Sigma_{i<k}\frac{p_i^T A r_k}{p_i^TAp_i}p_i \tag{5.32}$$  第$k+1$次迭代时取得的预测解向量计算公式为：$$x_{k+1}=x_k+\alpha_kp_k \tag {5.33}$$  其中$\alpha_k$的计算公式为：$$\alpha_k=\frac{p_k^T(b-Ax_k)}{p_k^TAp_k}=\frac{p_k^Tr_k}{p_k^TAp_k} \tag {5.34}$$  最终，共轭梯度法求解线性方程组的伪代码如下图3所示。图3 共轭梯度算法   然而，为了进一步加快求解的收敛速度，我们采用的预处理的共轭梯度法（Preconditioned conjugate gradient method，简称为PCG）$^{[2]}$，该算法预先定义了一个特殊的矩阵$M$，具体的算法流程如下图4所示。图4 预处理共轭梯度算法  预处理共轭梯度法等价于共轭梯度法计算如下线性方程组的解：$$E^{-1}A(E^{-1})^TE^Tx=E^{-1}b \tag {5.35}$$  其中，预处理共轭梯度法里面的矩阵$M=EE^T$，$M$也是一个保持不变的对称正定矩阵（上面的内容来自维基百科，我感觉这里是的矩阵$M$就是矩阵$A$）。目前已经有多种预处理方法，Robert Brdision在书中介绍了一种预处理方法，该预处理方法被称为不完全的Cholesky因子化（Incomplete Cholesky factorization）$^{[3]}$。  Cholesky因子化是这样的一个数学含义：给定一个正定矩阵$A$，可以找到一个这样的下三角矩阵$L$，使得$A=LL^T$。不完全Chkolesky因子化就是寻找一个近似于$L$的矩阵$K$，使得$A=LL^T\approx KK^T$。这里之所以不寻找确切的下三角矩阵$L$，是因为即便$A$是一个包含大量0元素的稀疏矩阵，但是其因子化的下三角矩阵$L$却不一定是稀疏矩阵，大规模的非稀疏矩阵将带来非常严重的内存问题。因此，我们希望寻找这样的一个因子化矩阵$K$，矩阵$A$中为0的元素，在矩阵$K$相应的地方亦为0，构成类似的一个稀疏矩阵。这就是不完全的Cholesky因子化。  上面讨论的不完全Cholesky因子化只允许矩阵$k$非零元素与矩阵$A$的非零元素具有一样的分布，这就是0级的不完全Cholesky因子化，通常记为IC(0)。首先我们将矩阵$A$分成一个下三角矩阵$F$和对角矩阵$D$：$$A=F+D+F^T \tag {5.36}$$  然后IC(0)的因子化矩阵$K$有如下的形式：$$ K=FE^{-1}+E \tag {5.37} $$  其中矩阵$E$是一个对角矩阵，因此对于因子化矩阵$K$，只有对角矩阵$E$未知，剩下的矩阵$F$可以直接从矩阵$A$中获取。Robert Brdison直接给出了计算对角矩阵$E$的计算公式，二维和三维的计算公式如下所示：$$E_{(i,j)}=\sqrt{A_{(i,j),(i,j)}-(A_{(i-1,j),(i,j)}/E_{(i-1,j)})^2-(A_{(i,j-1),(i,j)}/E_{(i,j-1)})^2}\\ E_{(i,j,k)}=\sqrt{A_{(i,j,k),(i,j,k)}-(A_{(i-1,j,k),(i,j,k)}/E_{(i-1,j,k)})^2 -(A_{(i,j-1,k),(i,j,k)}/E_{(i,j-1,k)})^2\\ -(A_{(i,j,k-1),(i,j,k)}/E_{(i,j,k-1)})^2} \tag {5.38}$$  然后Robert Bridson进一步指出可以在不完全Cholesky因子化的基础上进一步做一些提升，引入了Modified Incomplete Cholesky因子化，简称为MIC因子化。只需要在原来的IC因子化的基础上做一些简单修改，算法的收敛复杂度就从原来的$O(n)$降到了$O(n^{1/2})$，假设我们的欧拉网格宽度为$n$个格子。MIC因子化不再直接丢弃在A中相应位置为零的非零矩阵元素，而是考虑将这些非零元素加到对角线的元素上。  MIC因子化同样地构建一个下三角矩阵$K$，只不过此时$KK^T$每一行矩阵元素之和等于矩阵$A$的每一行元素之和。矩阵$K$同样是公式$(5.37)$的形式，只不过对角矩阵$E$的计算公式稍微有点差异：$$E_{(i,j)}=\sqrt{ A_{(i,j),(i,j)} -(A_{(i-1,j),(i,j)}/E_{(i-1,j)})^2 -(A_{(i,j-1),(i,j)}/E_{(i,j-1)})^2\\ -A_{(i-1,j),(i,j)}A_{(i-1,j),(i-1,j+1)}/E^2_{(i-1,j)} -A_{(i,j-1),(i,j)}A_{(i,j-1),(i+1,j-1)}/E^2_{(i,j-1)} } \\ E_{(i,j,k)}=\sqrt{ A_{(i,j,k),(i,j,k)}-(A_{(i-1,j,k),(i,j,k)}/E_{(i-1,j,k)})^2\\ -(A_{(i,j-1,k),(i,j,k)}/E_{(i,j-1,k)})^2 -(A_{(i,j,k-1),(i,j,k)}/E_{(i,j,k-1)})^2\\ -A_{(i-1,j,k),(i,j,k)} \times (A_{(i-1,j,k),(i-1,j+1,k)}+A_{(i-1,j,k),(i-1,j,k+1)}) /E^2_{(i-1,j,k)}\\ -A_{(i,j-1,k),(i,j,k)} \times (A_{(i,j-1,k),(i+1,j-1,k)}+A_{(i,j-1,k),(i,j-1,k+1)}) /E^2_{(i,j-1,k)}\\ -A_{(i,j,k-1),(i,j,k)} \times (A_{(i,j,k-1),(i+1,j,k-1)}+A_{(i,j,k-1),(i,j+1,k-1)}) /E^2_{(i,j,k-1)} } \tag {5.39}$$  MIC背后的数学原理是傅里叶分析，这里就不仔细展开了。在实际的实现中，为了有更佳的性能表现，我们混合使用IC因子化和MIC因子化，给MIC中增加的项乘上一个权重，通常取$0.97$或者更大。下图5展示计算对角矩阵$E$的伪代码，除了用e存储对角矩阵$E$之外，我们用precon存储$E$的开根倒数。图5 计算对角矩阵E   计算得到了precon之后，我们就将其应用预处理共轭梯度算法当中，这个预处理得到的下三角矩阵$K$被用于图4伪代码中的$z$向量，图4的伪代码中直接写出$z=M^{-1}r$，这是为了表示求解方程$Mz=r$。实际上我们并不是直接求$M$的逆矩阵，而是将$M$因子化成三角矩阵的乘积$KK^T$，变成了求解$K(K^Tz)=r$，因为$K$和$K^T$分别为下三角矩阵和上三角矩阵，因此我们可以分成两步来求解：首先求解$Kq=r$，直接前向替换即可；然后求解$K^Tp=q$，直接后向替换即可。具体的伪代码如下图6所示。图6 应用因子化的矩阵5、不可压缩投影  最后，流体的不可压缩投影$project(\Delta t, \vec u)$算法流程如下所示：Calculate the negative divergence $b$ （the right-hand side）with modifications at solid wall boundaries.Set the entries of $A$ (stored in Adiag, etc.).Construct the preconditioner (either MIC(0) for single threaded solves, or a more complicated domain decomposition set up for parallel solves).Solve $Ap=b$ with PCG（Preconditioned conjugate gradient）.Compute the new velocities $\vec u^{n+1} $according to the pressure-gradient update to $\vec u$.参考资料：$[1]$ Neumann boundary condition. From Wikipedia, the free encyclopedia$[2]$ Conjugate gradient method. From Wikipedia, the free encyclopedia$[3]$ Incomplete Cholesky factorization. From Wikipedia, the free encyclopedia$[4]$ 《Fluid Simulation For Computer Graphics》, Robert Bridson.

流体模拟概述 流体动画概述背景我目前致力于研究游戏级别的实时流体模拟，也可以说是流体动画。这个领域的知识内容比较多，我将把学习过程中收集到的资源和自己的心得记录一系列笔记，主要是流体仿真和渲染的技术，既可以方便自己回顾，也希望这些内容能帮助到大家。引言基于物理的计算机动画成为目前一个研究热点, 是计算机图形学发展的一个重要方向，而计算流体力学( Computational Fluid Dynamics, CFD) 也一直是热点话题。随着人们对电影和游戏的视觉效果越来越高，观众也希望看到真实流体的动画效果，传统的逐帧绘制的方式变得不再可行，于是基于物理的流体模拟变为主流。基于物理的方法都是通过求解运动方程进行模拟物体运动，而流体模拟中用到的运动方程就是Navier-Stokes方程。尽管该方程早在数百年前就已经提出，但它是一个非线性偏微分方程，手算结果几乎不可能，只能通过数值方法求解，随着计算机技术的发展，诞生了一门新的学科：计算流体力学，说简单点就是通过计算机庞大的计算力来解流体方程。流体动画的思想和方法就是来源于计算流体力学，但是两者的目标却不同，计算流体力学要求计算结果尽可能的精确符合现实，从而能指导工业设计。而流体动画则偏重于生成能够以假乱真的视觉效果，对精度要求并不高，但是计算速度一定要快，以达到实时播放的帧率。我侧重于实现流体动画，因此后面都是基于流体动画的实现展开。模拟方法传统的N-S方程求解方法通常分为两类: 基于网格模型的欧拉法和基于粒子模型的拉格朗日法。欧拉法从研究流体所占据的空间中各个固定点处的运动着手, 分析被运动流体所充满的空间中每一个固定点上的流体的速度、压强 、密度等参数随时间的变化,以及研究由某一空间点转到另一空间点时这些参数的变化, 该方法是一种基于网格的方法。拉格朗日法从分析流体各个微团的运动着手, 即研究流体中某一指定微团的速度 、压强、密度等描述流体运动的参数随时间的变化, 以及研究由一个流体微团转到其他流体微团时参数的变化, 以此来研究整个流体的运动, 该方法是一种基于粒子的方法。基于拉格朗日法的示意图：拉格朗日法中最常用就是SPH(Smoothed Particle Hydrodynamics)法了，它衍生出许多改进算法：总结本文先介绍一下级别概念，后面逐步细化每种算法。本文的信息可能还不够全面，会随着我学习的深入不断完善。

物理引擎基础 基础理论数学基础∇：向量微分算子、哈密尔顿算子、Nabla算子、劈形算子，倒三角算子是一个微分算子： $\nabla = {\frac{\partial }{\partial x}}\mathbf{i}+ {\frac{\partial }{\partial y}}\mathbf{j}+ {\frac{\partial }{\partial z}}\mathbf{k}$grad F=▽F，梯度（gradient），梯度指向标量场增长最快的方向，梯度的长度是这个最大的变化率： $▽f=\frac{\partial f}{\partial x} \mathbf{i} + \frac{\partial f}{\partial y} \mathbf{j} + \frac{\partial f}{\partial z} \mathbf{k}$div F=▽·F,散度(divergence),描述了通量源的密度，可用表征空间各点矢量场发散的强弱程度： $\operatorname{div}\,\mathbf{F} = \nabla\cdot\mathbf{F} =\frac{\partial U}{\partial x} +\frac{\partial V}{\partial y} +\frac{\partial W}{\partial z }$Δ、$∇^2$ or ∇·∇：拉普拉斯算子（Laplace operator），定义为梯度（▽f）的散度（▽·▽f）： $\Delta f = \nabla^2 f = \nabla \cdot \nabla f$, $\Delta f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2}$rot F 或curl F=∇ × F,旋度(curl，rotation)，是向量场沿法向量的平均旋转强度: $\operatorname{curl}\,\mathbf{F} = \nabla \times \mathbf{F} = \left(\frac{\partial F_3}{\partial y}- \frac{\partial F_2}{\partial z}\right)\mathbf{e}_1 - \left(\frac{\partial F_3}{\partial x}- \frac{\partial F_1}{\partial z}\right)\mathbf{e}_2 + \left(\frac{\partial F_2}{\partial x}- \frac{\partial F_1}{\partial y}\right)\mathbf{e}_3$基本关系：一个标量场f的梯度场是无旋场，也就是说它的旋度处处为零：$\nabla\times (\nabla f) = 0$一个矢量场F的旋度场是无源场，也就是说它的散度处处为零：$\nabla\cdot (\nabla \times \mathbf{F} ) = 0$F的旋度场的旋度场是： $\nabla \times (\nabla \times \mathbf{F} ) = \nabla(\nabla\cdot \mathbf{F}) - \nabla^2 \mathbf{F}$

收集的物理引擎资源 收集的物理引擎资源YangWC's Blog飞翔的子明游戏物理笔记何小伟的网站Ten Minute PhysicsVideo Game Physics Tutorialwmdcstdio的博客CFD资料汇总未完待续~