流体仿真-基础篇

2018-11-14 | 阅读：次

前言

主要是根据Robert Bridson著名的Course Notes，个人写的笔记整理，补充上一些我认为不够清晰的地方。可能分成几个部分，首先先整理基础篇，从基本的流体力学方程入手，解释每一项的具体含义，以及原理，然后针对不同的流体特性对方程进行简化，而且NS方程目前还没有人能够求出解析解，只有100多个特解，因此只能采用数值上的近似方法，所以接下来会进入数值分析该方程的步骤，近似求解方程，解出速度场，然后进行流体的模拟仿真。需要有一定的数学基础：微积分，向量分析，微分方程数值分析等。

流体力学方程

要进行流体仿真，首先就要搞清楚流体所遵循的物理定律，也就是力学方程，然后才能进行准确的模拟。流体力学中关于流体最重要的方程就是纳维-斯托克斯方程，也就是NS方程，具体定义如下：

该方程由两部分组成，一式是动量方程，二式是质量守恒方程，接下来一一解释。在解释之前，先把符号含义定义，u向量指的是流体的速度，p是压强，g向量是重力加速度， ρ 是流体的密度，希腊字母ν称为运动粘度。

   由于表达不方便，文中的大部分向量的箭头标志我都省略，比如u，a，g，F等等。

动量方程

动量方程指的就是一式，乍一看很复杂，其实就是牛顿第二定律，也就是外力的大小等于动量的变化率，或者说F=ma。首先a就是加速度，应该等于速度对时间的导数Du/Dt，这里的大D指的是全导数，在流体力学里，这个导数也叫作物质导数或者随体导数。后面会说到这个随体导数其实是在拉格朗日方法中用到的，随体就表示跟随物体，先不讨论。

根据牛顿第二定律，ma=F，a用随体导数替换，也就是mDu/Dt=F，m是流体微团的质量，这里用微团表示我们所关注的流体小部分，这是流体力学中经常用到的术语。F就是这个微团所受的外力，那么显然，这个外力具体是什么就是接下来的目标。

首先比较清晰的是重力，也就是mg，这个微团肯定是会受到重力的
然后是压力，不同的流体微团之间一定会存在压力，压力高的区域挤压压力低的区域，解决压力的问题也是流体仿真中非常非常重要的一部分。那么压力是什么呢？压强其实就是压力场的差，压力场的差也就是压力场的梯度，所以我们可以简单的用压力场的梯度的负方向−∇p表示压强，压力就是压强在这个流体微团体积上的积分，一般为了简单起见，而且我们假设流体微团足够小，可以直接乘以这个体积近似整个微团的压力，也就是-V*∇p
然后要考虑的就是粘性力，流体之间是有粘性的，区别只是粘性的大小不同，比如蜂蜜的粘性很大，平时的自来水粘性小一些。为了数学建模，粘性从另一方面看其实就是为了让当前的流体微团和周围的流体微团的速度相似，即周围流体微团的平均速度。这用拉普拉斯算子来描述在准确不过了，拉普拉斯就是描述了与周围平均速度差别的算子，为0时，和周围平均速度一致。拉普拉斯算子记为∇ · ∇，也可以记为∇的平方，其实也就是梯度的散度，展开来说就是空间的二阶导数。在向量分析领域，拉普拉斯算子有很重要的意义。和压力一样，为了得到整个微团的粘度，必须对体积积分，也可以简单的乘以体积近似，即V*∇ · ∇u。但是这样衡量不了不同粘性流体的区别，因此还需要加上一个动力粘性系数 µ。

把上面的三个外力组合起来可以得到整个流体微团的动量方程：

这个方程要想没有错误，必须取极限，也就是微团数目趋向于无穷多，并且每个流体微团的大小尺寸趋向于0才行，不然这个方程就是错误的。现在有一个问题出现了，如果取了极限，那么为微团的重量m和体积V也同样趋向于0了，这就无法计算了，我么可以方程两侧同时除以V，把m和V替换成密度：

现在看起来和一开始给出的NS方程中的动量方程就有点相似了，左右再同时除以密度并且整理一下得到：

为了进一步简化方程，定义ν = µ/ρ，将动力粘性系数除以密度的商定义为动力学粘度ν，得到：

现在看来和NS动量方程更加接近了，唯一的区别就是Du/Dt，也就是我们先前说到的随体导数，下面我们将理解随体导数具体是什么，有什么意义。

拉格朗日描述和欧拉描述

在前面说到随体导数的时候我说道这是在拉格朗日法中用到的，指的就是这里的拉格朗日描述。何为拉格朗日描述？就是把流体看成一个个的粒子组成，也可以吧前面的流体微团看成粒子，每个粒子单独拥有自己的属性，包括位置x，速度u等，这个描述和我们对流体的直观理解更加接近，实现起来也较为容易，但是缺点是流体粒子之间的交互变得复杂了起来，因为每个粒子都是独立的。而随体导数描述的就是一个量随着时间的变化率，这里不特地指定这个量是速度u，更加一般化了，也可以是温度之类的，恰好符合拉格朗日的描述，而拉格朗日法的一大好处就是由于粒子的数目是固定的，因此NS方程中的质量守恒方程是自动满足的。

欧拉描述则是一个新的视角，他不把流体看成粒子系统，不会追踪每个粒子，而是只观察空间中的特定的固定点，观察这个点的属性随时间的变化，而这个固定点属性的变化是由流体流过而改变的，欧拉描述和人们的直观印象相对不是很符合，但是他也有很多的优点：

更容易处理空间导数，比如压力梯度和粘性的拉普拉斯算子，只需要利用空间点的有限差分即可，而粒子系统比较混乱，难以计算
欧拉网格上近似空间导数比随机移动的粒子方便多了

因此不同于拉格朗日描述遵循的质量守恒，欧拉描述遵循的是体积守恒。

其实还有一种LBM法，从微观力学的角度考虑，但是笔记中没有提到，也不再多说

而拉格朗日法和欧拉法就是通过随体导数相关联的。为了一般化，我们取一个一般的量q，他关于时间t和位置x，即q(t,x)，不同时间不同位置q值不同。现在求q对时间的全导数：

我们可以发现，q对时间点额全导数就是随体导数，这在开始我也提到了，他是在拉格朗日描述中用到的，展开后发现包含了q对时间的偏导，这是与x无关的，也就是说x是常数，这个偏导代表了q和常数x位置随时间的变化，正好满足了欧拉描述。而剩下一项是一个修正项，与流过这个固定点的流体有关。组合起来就是拉格朗日描述就是欧拉描述加上经过欧拉网格点的流体带来的属性变化。

为了完整，我给出随体导数在三维上的完整形式：

随体导数可以用来表示平流或者说对流，在流体力学中随体导数的那一项也叫作对流项。所谓对流，或者说平流，就是一个量在一个速度场中被速度场输送的过程。这一般很好理解，但是当这个量也是速度的时候可能会产生一些混乱，不仅是温度之类的可以再速度场中平流，速度本身也是可以的。一般我们取随体导数为0，就是一个简单的平流方程，他表示这个量仅仅在速度场中到处移动，在拉格朗日的视角下，他没有任何变化，即他的温度，速度等没有变化。但是在欧拉视角下，可能发生了变化。可以举一个简单的例子：

  一个一维空间上，不同位置的温度不同，满足T(x)=10x,即x=10处的温度为100.现在这个空间置身于一个流速稳定的速度场中，也就是u=c，是一个常量。
  令DT/Dt=0，解这个平流方程(全导数展开即可)，得到T对t的偏导数为-10c，这意味着对于欧拉描述来说，空间上的某个固定点的温度和速度场的速度大小和
  方向有关。如果c=0，则每一处的温度不变，若风从正方向以c=2吹过，这点的温度将以20的速度上升，下降亦然。

不可压缩性

自然界中可压缩流体的仿真很复杂，而且很不常见，所以我们可以只考虑不可压缩流体。而且即使是可压缩的流体，也不太可能发生很大体积的改变，因此我们可以视为不可压缩。不可压缩，另一个角度来看就是体积不会改变，密度不会改变，假设一个流体微团，表面为∂Ω，体积为Ω，则我们可以定义这个微团的体积变化率为：

体积的变化率等于每小块表面速度法向部分在整个面上的积分，面积分结果大于0，代表总体体积变大，向外的速度大于向内的速度，微团被拉伸。为了得到不可压缩流体，因此需要满足这个面积分的结果为0.

再根据散度定理，把微团面积分转化为整个微团的体积分(关于散度定理不展开，就是一种将面积分转化为体积分的方式，类比格林公式)：

观察到对任意的Ω而言，也就是流体的任意区域，这个方程都应该成立，唯一满足这个条件的函数就只能是0了，因此得到∇ · u = 0，注意到这就是NS方程中的二式，我们一开始称为质量守恒方程，其实同时也是不可压缩条件，正因为不可压缩，所以体积不变，然后才能保证质量也不变。满足这个不可压缩条件的向量场我们称为“散度自由”或者“零散度”。事实上，仿真不可压缩流体最重要的就是要保证流体时刻处于零散度状态。零散度也是压力产生的原因，因为体积保持不变，因此微团一面受力，这面体积缩小，为了满足体积不变，另一面体积就要增大，这就产生了对另一面的压力。

舍去粘度

实际上对于大部分场景，粘度项是可以舍去的，一是因为常见的场景下粘度是可以忽略的，现实大部分流体都可以看成非粘性不可压缩流体；二是因为我们会用数值的方法来求解NS方程，这实际上会引入叫“数值耗散”的问题，也就是误差，后面也会分析这个误差其实就是关于速度的空间二阶导数，空间二阶导数其实就是拉普拉斯算子，所以实际上不可抗的我们都会引入粘度项，或者叫类粘度项，因此公式中的粘度可以显式的舍去。

舍去粘度项的NS方程，我们叫做欧拉方程：

接下来我们主要解决欧拉方程，使用欧拉描述的方式，后面也会介绍SPH这种拉格朗日描述的算法。

边界条件

流体仿真最重要的一部分是边界条件，在计算流体力学这个领域中，大部分的精力都花在边界条件上。我们之前一直在讨论流体的内部情况，没有涉及流体的边界问题，比如流体和固体墙的交互，再比如流体表面和空气的交互，这两个也是最重要的交互，还有一种是不同流体之间的交互，比如蜂蜜和水的交互，这个就复杂很多了，暂时不在考虑范围内。

所以我们主要讨论“固体墙”和“自由面”两种情况：

固体墙

流体和固体墙的交互，有两种情况，一种固体墙是完全固定的，那么为了不让流体穿过固体，可以简单的将边界的速度法向分量设置为0，即：

另一种情况是固体是移动的，也就是固体本身也有一个速度，那么我们就可以令边界法向速度分量等于固体的速度，让流体跟随固体而走，我们没有限制切向速度，因此流体可以自由的沿着固体边界切向滑动，这也是符合实际情况的。

 注意到我们讨论的非粘性流体，如果是粘性流体的话，固体墙的粘性对流体的切向的速度分量也是有影响的

另外对于某些特殊的固体墙，比如通风口或者排水管，我们需要特殊处理，不会直接让速度法相分量等于固体的速度，而且赋予抽水泵汲取或者喷出的速度。

#### 自由面

自由面就是我们终止流体模拟的地方，比如水的表面，这时候需要考虑空气，但是因为空气比水轻700多倍，因此可以把空气看成拥有大气压强的一片区域，在实际仿真中我们只关心压力差，可以简单的把大气压设置为0，即自由面是压强为0的面。

液体表面有时还有一种特殊情况需要考虑，就是水的表面张力。从分子层面来看，表面张力的成因是因为液体分子之间的吸引力大于水分子和空气各分子之间的吸引力造成的。为了数学上的建模，可以从几何角度来看，表面张力就是尽力最小化水的表面面积，或者换一种说法，尽量减少表面的平均曲率，我们使用第二种表述来进行建模。

如下公式是成立的，表面之间的压力跃变正比于表面的平均曲率，即：

其中的 γ 是表面张力系数，已经有精确的物理测量结果了，常温下，水和空气的γ ≈ 0.073N/m。换种说法， you表面张力的自由面，书表面的压力就是空气压力(通常假设为0)加上压力跃变。

数值近似

上面已经基本很详细的剖析了NS方程各项的意义和背后的数学建模原理，以及针对常见流体，也就是非粘性不可压缩流体的简化形式，然而即使经过简化，变成了欧拉方程，解决起来依然不是很简单，无法求出解析解，因此流体仿真一般采用数值近似的方式。

数值近似直接解决整个欧拉方程是很复杂的，因此一般采取分治法的思想，分而治之，把欧拉方程分离成几个部分分别解决，然后再结合起来，这在数学上是合理的，假设：

这里f和g函数可以是任意的黑盒子函数代表分离的部分，我们可以将这个微分方程使用显式欧拉分离成：

这个分离是一个一阶精确算法，可以泰勒展开验证一下：

正如你看到的这样，原来的微分方程被分离成两个单独的方程，一个只包含f函数，一个只包含g函数。需要注意的是分离后的方程的解决顺序也有讲究，一般顺序解决，因为有的积分器需要另一个积分器的输入，顺序序列能保证不出错，并行计算可能发生问题。

 显式欧拉：也叫前向欧拉，是一种常用的微分方程数值近似方法，利用下一步的函数值等于当前函数值加上
 时间步delta乘以当前函数的导数这个关系，f(x+delta) = f(x) + f'(x)*delta,因为向前看了一步，
 所以叫前向欧拉，自然相应的还有后向欧拉，就是向后看一步，后向欧拉
 是前向欧拉的改进方法，是稳定的算法，而前向欧拉是有不稳定区域的，但是后向欧拉复杂了很多。

分离欧拉方程

我们同样可以分离欧拉方程，分别解决，一般分为平流部分，体积力部分和压力/不可压缩部分：

分开解析一下：

平流部分，可以表示成advect(u, ∆t, q)，表示一个量q在速度场u中输送一个∆t的时间步，就是我们之前说到很多次的平流方程。
体积力部分，此处只考虑有重力的情况，浮力之类的体积力也可以合成到重力中。公式表示了重力加速度对流体速度的影响，可以用前面解释的显式欧拉法来近似，这是最简单的一部分。
压力部分，表示成project(∆t, u) ，将压力梯度移到方程的右侧就可以看出，这是考虑了压强对速度的影响，同时还要保证满足不可压缩条件，这也是保证压强的来源。

一旦我们上述分离后的方程都解决了之后，只需要按如下方式组合起来即可:

在进入具体方程的数值近似之前，我们先考虑几个问题：

时间步的选择，数值近似就是离散化的过程，时间步决定了离散的精度，因此时间步的选择是一个trade-off问题，要兼顾速度和准确度。有特殊需要注意的地方，就是时间步不能越过当前动画帧的间隔，也就是不能tn + ∆t > tframe，一般超过，我们需要强制令∆t = tframe − tn，同时设置一个标志标明我们到了这一帧的结尾，这是为了进行一些操作，因为每一帧的结尾我们需要渲染到屏幕或者是保存状态。
空间离散。之前包括∆t，前向欧拉等等都是在时间上的离散，一个个时间步，现在我们考虑空间上的离散，因为我们使用的是欧拉方法，因此空间离散就使用网格(Grid)，以前的很多网格法都是规则网格，也就是所有的数据都存在网格的中心，包括速度，压力等等。这种网格会有一些很大的困难，比如解决压力梯度时，会出现null-space，也就是零空间的问题，这个后面会说。我们采用的是MAC网格，这是一种交错网格，数据交错存储，不是全存在网格的中心，下面给出两张图，分别是2D和3D的情况，这两张图需要牢记：

可以看到不同的变量保存在不同位置，所以称为交错网格。这种结构的好处就是方便利用差分计算压力梯度和速度场的散度，接下来在相应的计算部分会详细介绍。

 差分计算：有必要先介绍一下差分计算，这也是常见的一种数值近似的方法，和前面的前向欧拉法一样，不过前向欧拉一般用于时间上的近似，差分适合空间上的近似，
 差分法又分为前向差分，后向差分和中心差分。举个例子，一维上三个点pi-1,pi,pi+1，求pi出的导数，可以用差分近似，前向差分：(pi+1 - pi)/delta ，后向差
 分：(pi - pi-1)/delta ，中心差分：(pi+1 - pi-1)/2\*delta，其中前向后向差分有偏差，精确到O(∆x)，中心差分无偏差， 而且精确到O(∆x2)，所以一般使用
 中心差分。
 
 但是中心差分有一个问题，就是完全不考虑pi本身的值，中心差分计算为0，本来应该是导数为0，也就是说函数是常数，但是中心差分不考虑pi，因此pi可能在任意位
 置，因此不能保证是常数函数，这就是前面说到的null-space问题，零空间问题，这常常是致命的。
 
 解决方案就是使用MAC网格，注意到MAC网格使用中心差分(pi+1/2 - pi−1/2)/2\*delta，因此没有跳过任何值，所以没有零空间问题。
 
 关于高阶差分问题，这里也简单提一下，一阶中心差分是：p'(i) = (pi+1 - pi−1)/2\*delta,二阶差分公式是：p''(i) = (pi+1 + pi−1 - 2\*pi)/delta平方,至
 于是如何得到的，主要是利用泰勒展开，具体是在x+deltax和x-deltax处泰勒展开，然后俩式相加整理一下就可以得到二阶差分公式，更高阶的也是如此。

MAC网格主要的缺点是求解空间某处的完整的速度时一定需要插值计算，因为MAC网格把速度各个分量存在了不同的面上，参考上面的图。

平流方程

接下来解决实现具体方程的数值近似了，体积力的近似简单的用前向欧拉近似一下就可以了，重点是平流方程和压力方程，这两个解决了，欧拉方程就解决了，欧拉方程解决了，流体仿真问题就基本解决了。

平流方程就是随体导数为0，前面也说过很多次，表示成advect(u, ∆t, q)，表示了给定一个速度场u(经过MAC网格离散过得)和δt和当前的场量qn，求解下一时刻q在速度场中平流输送后的结果。

首先考虑的就是直接暴力解微分方程。

直接解微分方程，然后有限差分近似导数，比如：使用显式欧拉计算时间偏导，然后利用中心差分计算空间偏导，然后上式可以近似转化为：

整理一下可以得到：

我们看起来好像获取了场量q的更新方式，但是实际上有问题的。前向欧拉在这里是无条件不稳定的，我们之前也提到过，显式欧拉是不稳定的，隐式欧拉是稳定的，在这里如果我们使用显式欧拉的话，不管∆t取多小，都会发生数值爆炸。数学上是因为中心差分法生成的雅可比行列式的特征值是纯虚数，始终在欧拉方法的稳定区域之外。而且即使解决了这个问题，空间导数也存在之前说过的零空间(null-space)问题。

所以暴力解决微分方程是有困难的，但是也不用担心，很多的数值方法对这些问题都有很好的分析，而且提供了很多的针对空间导数的更复杂更有效的有限差分方式。

但是我们重点将放在另一种完全不同的方法上，叫做半拉格朗日法，这种方法更简单，更容易实现，而且更基于物理的动机。我们注意到如果使用拉格朗日描述的话，平流方程其实完全不用担心，他是可以被自动解决的，尤其是我们使用粒子系统的话。为什么呢？因为我们的目的是解决场量q在平流输送后下一时刻的更新值，而在拉格朗日描述的粒子系统中，这个值就是此刻经过这个位置的粒子在上一个位置时的旧值。应该很好理解，这个粒子流动到这个位置带来了新的值，因此我们只需要解决两个问题：

如何找到上一个位置？
上一个位置的值是什么？

找上一个位置，其实就是逆向速度场，将当前的粒子逆向移动一个时间步，然以后计算这个位置就行了。举一个例子; 想象一个粒子，从位置xp开始，遵循常量速度u在速度场中移动，经过一个时间步∆t停止在xg点，那么现在我们逆向速度场计算一下，xp = xg-∆t*u，可以看到就是简单的使用了一步显式欧拉，显示欧拉法一般的确足够了，但是改进的欧拉法可以取得更好的效果，比如一种改进欧拉法，基于二阶龙格库塔(RK2)，这里不多说。

经过上面的例子，我们得到了上一步的位置，下面自然就是解决那个位置具体的q值了。我们知道，由于我们使用的是欧拉法，具体也就是MAC网格，那么我们回退得到的位置其实很大程度上都不会在网格点上的，因此我们一般使用线性插值来通过周围网格点上的q值来进行计算，一般2D使用2次线性插值，3D使用3次线性插值：

注意这个q值代表的是一般的场量，可以具体是压强，烟雾的密度，温度，速度等等，这些不同的变量中，一般只有速度时存储在交错网格未知的，也就是各个面上，其他
变量一般都是存储在中心点的，可以回忆之前的那两张MAC图。

PS：这个粒子实际上是不存在的，是我们虚构的，我们只是利用这个一个虚构的粒子来直观的解决平流问题，也正是因为我们在欧拉法中使用了一种类似粒子系统的拉格
朗日思想方式，所以这个方法才被叫做半拉格朗日法。

边界问题

需要注意一个问题，就是我们使用半拉格朗日法的时候，回退的粒子可能回退到比如水的外面，回退到空气中(自由面的外面) 或者固体中，那么我们进行插值时就不能简单的取周围网格点的值进行了，我们就需要考虑边界条件。注意回退到固体中一般是因为水是从那里流入的。

水是从固体中流出的，那么当我们会退到固体中的时候，这个值就不能用插值进行计算，但是有个条件就是我们必须知道流入的具体的量，比如：水从一个区域流入，以特定的速度v，那么半拉格朗日法回退到这里时就令q等于v，不需要插值，而且这里的量比如v，我们一定要已知。(当然也可以是其他的量）
水回退到自由面的外面，可以看成是回退到水面以外，那么自然也不能用周围网格进行插值，因为周围网格包含空气网格。这时候我们需要用周围在自由面上的邻近的网格点外推这个量。

上述两个边界问题，一个是由于流体从外部流入引起，另一个可能是由于误差引起，这个误差可能是显式欧拉或者RK2算法计算位置时带来的计算误差。

时间步

半拉格朗日法回退粒子需要定义一个回退的时间步，那么时间步的选择就需要考虑了。在数值模拟中，稳定性是一个重要的考量，不会发生数值爆炸很重要，因此需要谨慎选择时间步，但是对于半拉格朗日法而言，这个问题不需要担心，因为拉格朗日法的回退的值都是通过插值计算的(边界条件特殊处理)，因此不会出现数值爆炸的，一定是有一个边界的，因此我们可以把注意力放在准确性和速度的衡量上，不怎么需要担心数值爆炸的问题。

当我们令时间步等于帧长的时候，速度很快，但是准确性不行，相反，步长取小，准确性高了，速度就很慢了，因此这也是一个trade-off，幸好很多先辈已经给出了时间步长的建议区域，一篇论文给出了：

其中∆x是网格的边长，umax是流体速度最大值的一个估计，关于umax更加鲁棒的估计还考虑了体积力的影响，比如重力或者浮力：

这个改进还有一个好处是保证永远为正，避免了除零错误。

数值耗散

终于到了解答前面一个疑问的时候了，前面在丢弃粘性力的部分提到过，即使直接丢弃了粘力，由于数值耗散的原因，其实是迫不得已已经考虑了粘性力。这里给出详细的解释。

首先说明一下数值耗散产生的原因，我们在半拉格朗日法中，计算旧值需要进行线性插值，我们使用的是前一个时间步周围值的加权平均来计算的，所以每一步平流，我们都做了一次平均操作。而平均操作会平滑尖锐的特征，这我们就称之为耗散。

用一个具体的例子来说明，为了简单，我们取一个一维的平流问题：

这个一维的速度场是一个常数u，并且u>0，我们假设∆t < ∆x/u，这代表我们回退一步粒子的长度不会超过∆x，也就是不会超过这个网格单元。那么假设我们现在的位置是xi，前一个位置是xi-1，那么我们在xi处回退粒子的时候，粒子新的位置还是[xi−1, xi]，因为最近的网格点就是xi和xi-1，因此我们线性插值计算一下新的q值：

整理一下：

观察一下这个式子，其实就是用显式欧拉计算时间导数，用前向差分计算空间导数，我们把qi-1泰勒展开：

然后代入上述方程：

同样观察这个式子，发现把qin和∆tu q对x的偏导部分移到公式左边，方程左右在同除时间步∆t，在二阶截断误差下，就得到：

公式左边就是我们的平流方程的主体，右侧不等于0，而是等于一个空间二阶导数，由此可见这个空间二阶导数就是带来的误差项，前面我们也说过了，空间二阶导数其实就是拉普拉斯算子的数学本质，梯度的散度就是空间二阶导数，所以说这个误差项就是一个一维的拉普拉斯算子，一个类粘度项。这就是我们为什么明明舍去了粘度项，实际上还好像有粘度的原因。这个现象就叫做数值耗散。可以考虑使用更优的插值方法来减轻数值耗散，但是一般代价很大，而且效果一般。

压力/不可压缩性

平流部分大体告一段落，接下来进入压力/不可压缩部分的数值近似，把这个解决了，欧拉方程的数值近似就解决了。压力部分包含两个部分，压力梯度部分和不可压缩性条件：

上述方程就是欧拉方程分离后压力部分的显式欧拉展开。

此外还要满足固体墙边界条件：

因此project(∆t, u）的主要任务就是从速度场u中减去压力梯度，并且同时需要满足不可压缩条件和固体边界条件。

因此我们的任务也就相应分成两部分：

离散化压力梯度(MAC网格中)
离散化MAC网格中的散度

我们暂时考虑的边界都是和网格对齐的边界，一个网格是空气或者是固体，就代表着整个网格里都是空气或者固体，当然这在实际模拟中肯定是不真实的，
因此后面会介绍不规则边界和部分填充网格的解决方案。

离散化压力梯度

压力梯度的离散化是很简单的，之前介绍MAC网格的优点时就说过，解决压力梯度很方便，我们把每个网格的压力都存储在网格单元的中心，速度存在各个面上，因此离散化压力梯度简单的使用有限差分就可以了，二维中的形式是:

三维上同理，把z轴包含进去就行了。我们可以看到压力在每个速度分量上进行更新，同时也可以看到公式依赖于外部压力条件，因此具体边界条件很重要，我们暂时只认为网格单元都是纯的，比如没有部分水部分空气的情况。

对于自由缅边界条件，就很简单，假设外部压力位0，也就是把公式中相应的p赋值为0即可，这种直接指定边界处值的边界条件我们一般称为狄利克雷边界条件。

还有一种固体墙的边界条件，这个就稍微复杂一些，这种被称为诺依曼边界条件。此时，边界处减去压力梯度的时候应该满足固体边界条件，也就是边界速度的法向分量应该等于固体速度的法向分量，也就是边界处下一时刻的速度等于固体的速度，所以可以用固体速度对之进行替换，原来的压力方程是：

替换后在重新整理一下，变成压力更新公式：

这个压力更新公式十分重要，这为我们提供了计算诺依曼边界压力的依据。这里可能会有一些疑问，就是我们计算固体边界的压力时，是通过这个压力更新公式计算的，那么这时候如果这个固体单元的另一面也和流体相邻，那么压力也是用这个公式的，问题就来了，两个压力不一样怎么办？因为我们MAC网格的压力是保存在中心的，因此压力值只有一个。这里就需要明确一点，对于流体单元的压力确实是这样的，但是对于固体单元的压力，我们实际上是不存储的，其中一个原因就是存储会产生冲突，另一个原因是没必要存储，因为我们可以通过压力更新公式计算出来，我们可以想象成固体单元的压力也是存在各个面上的，但是只是想象，数学上是虚构的。

你或许会想，这里解决边界速度在压力影响下的更新时，为什么不直接把边界速度赋予固体速度，还要先去算压力梯度，然后再减去这个压力梯度呢？

事实上，有很多方法确实是这样做的，但是之后依然还要来考虑压力梯度，因为当我们解决流体内部压力问题时，我们必须知道这个墙边界的压力是什么。

离散化散度

第二个任务就是离散化散度，散度公式是∇ · u=0，我们可以直接展开解这个微分方程，二维展开是(3D同理，可以自己展开看看)：

有限差分进行离散;

注意散度我们只估计流体的内部，固体不管

还有一种离散化的方式，估计流体对某个网格单元的流入率和流出率，这是从散度的物理意义上考虑的，如果流入和流出一样，那么这个单元的散度自然就是0。

在精确连续介质中，就是围绕网格单元面速度的法向部分的积分：

MAC网格前面也说到估计散度很方便，就是针对这种思路。因为MAC的速度以分量的形式存在各个面上，所以只需要简单的用法向速度乘以面的面积即可。尺度变换后你会发现和中心差分公式计算的结果是一样的。这种我们直接估计一个量在网格面周围积分，而不是关注微分方程的数值方法叫做有限体积法。

有限体积法具体过程：面积是∆x的平方，所以散度就是∆x的平方乘以格格面的速度分量之和，注意速度我们都是按照朝外为正，朝内为负。计算后你会发现，左右同时除
以∆x的三次方的话，右侧和直接有限差分散度公式的形式是一样的，左侧则尺度缩小了∆x的三次方。

这里如果我们不用MAC网格的话，会发现出现了null-sapce的问题。

压力方程

把离散化后的压力更新公式和散度公式结合起来就是最后的压力方程了，我们最终的目的是保证流体内部无散度，即：

上式就是散度公式的差分近似形式，也是上面离散化散度部分推导出来的，然后再把压力更新公式带入，得到：

三维也是同理，可以自行推导。观察这个方程，这就是压力方程的最终版，左边括号里的部分，之前在高阶差分近似的时候介绍过了，应该很熟悉，这就是有限差分二阶导数的形式的相反数，是由泰勒展开推导出来的，而二阶空间导数就是阿拉普拉斯算子，已经提到很多次了，而右侧括号里也就是散度的数值近似，因此这个公式其实就是泊松问题的数值近似。

拉普拉斯算子等于0的方程叫做拉普拉斯方程，也叫调和方程，它的解就是调和函数，在各个领域都有重要意义。
而右侧不为0的方程叫做泊松方程，拉普拉斯算子带系数的问题也叫泊松问题。

上述的压力公式只是一般形式，需要根据不同的边界条件，狄利克雷边界条件和诺依曼边界条件进行调整。需要把公式中相应的压力换成符合边界条件的形式。空气网格单元就把对应的p设置为0，固体网格单元的p利用前面的压力更新公式进行计算，再带入。

具体对于代码而言，我们一般做以下几件事：
1）对于气体边界条件，简单的从方程中删除那个p
2）对于固体边界条件，我们也删掉那个p，但是也降低pij前面的系数，pij前面的系数等于周围非固体网格单元的数量。
3）改变右侧测量的散度，使用固体墙的速度来代替旧的流体速度，ui+1/2，j。

矩阵向量形式

注意上面的压力方程，我们只是计算了第(i，j)个网格单元的压力情况，而一般而言网格单元可以成百上千，因此采用矩阵向量形式可以极大地提高计算效率。注意到压力方程左侧为压力相关，右侧为散度，可以用如下形式表示:

其中p就是压力向量，多少个网格单元就有多少个列，一一对应，d为散度向量，和p的维度一样，表示各个网格单元的散度，A是系数矩阵，每一行对应一个压力方程(因为每一行代表一个网格单元，每个网格单元可以计算一个压力方程)，仔细一想就知道这个系数矩阵A是很稀疏的，大部分的位置就是0，由于每一行代表一个网格单元，因此每一行对于2D来说其实最多只有5个非0值，就是压力方程左侧表示的上下左右四个和自己本身的压力系数的，同时也可以发现系数都是的倍数，因此可以方程左右同时除以这个数，那么系数矩阵就简化为非对角项部分只有0或者-1，对角项部分是 p(i,j)前面的系数，具体来说就是(i,j)周围非固体的网格单元数目。

这个A还有这更深的数学性质，比如说它是对称正定的，因此在二维上，每一个网格单元我们只需要存储三个数字，之前提到了每一行A最多在二维上只有5个非0值，因为对称，所以只需要存3个，即对角线元素，以及两个正方向上的相邻单元记录。我们一般记为：

二维上：Adiag(i,j), Aplusi(i,j) ，Aplusj(i,j)
三维上： Adiag(i,j,k), Aplusi(i,j,k), Aplusj(i,j,k) ， Aplusk(i,j,k)

共轭梯度法

解决想Ax=y这种正定线性方程系统，有很多经典的方法，共轭梯度法就是其中一种(CG)。共轭梯度法是一种迭代方法，优点就是操作简单，计算简单，写成代码也很方便，所以很受欢迎，但是缺点也就是网格越大，收敛越慢，因此我们会使用CG算法的改进，预处理器共轭梯度算法(PCG), 接下来一一介绍。

关于CG法的细节不多说，注意一般系数矩阵A为单位矩阵时迭代是最快的，比如Ip=d，那么解自然就是p=d，所以我们就想能不能把A改进成单位矩阵或者接近单位矩阵。所以就把方程左右同时乘以一个矩阵M，编程MAp=Md，如果M接近A的逆矩阵，那么MA就接近单位矩阵，这样一来共轭梯度算法的迭代就会很快。这种加了一个预处理器M的方法就叫做预处理器共轭梯度算法(PCG)。

如何判断收敛呢？

PCG算法解决Ap=d，令残差r=d-Ap，也就是散度减去散度，就是散度差，因此残差就是更新压力的估计之后，速度场剩下的散度值。如果这个r变成0了，
就代表剩下散度为0了，压力方程就得解了。

但是在实际操作当中，我们不需要让残差r正好为0，而是等于一个较小的数就可以了。另外为了防止不精确浮点计算导致无法收敛，一般还需要制定一个
最大迭代次数。

还有就是关于压力的初始化选择，一般实际操作中初始化为0就行了。

所以PCG具体算法流程为：

之前的分析中，令M=A的逆固然好，但是计算代价很大，所以我们一般选择其他的预处理器。这就要引入乔里斯基分解了。

不完全乔里斯基分解

乔里斯基分解可以看成对于一个正定对称矩阵而言，可以分解成一个上三角矩阵和下三角矩阵的点乘，即：

因此Ap=d就可以用L(LT p) = d来替换，然后我们可以分两步计算，先把LT p看成整体，计算Lq=d，然后再解决q，也就是LTp=q，其中T代表转置。

这种方法也带来了一些问题，那就是A是稀疏的，存储起来占用内存很少，但是L不是，因此会出现内存溢出。所以引入了不完全乔里斯基，这个方法解决了内存问题。

在算法要创建一个新的非零值，但是这个位置在A中是0的时候，把这个非零值设置为0.这样就没有内存问题了。但是问题是A和LLT不在相等，但是好消息是他们足够接近，因此可以近似使用。

假设我们把A分解成A = F + D + F T，其中F为下三角矩阵，D是一个对角矩阵，其实也就是把A中的各个项提取出来。那么因为A=LLT，因此可以解出L = FE−1 + E，-1代表逆矩阵。，E是对角矩阵。现在我们可以看出了，我们其实只要存储 L矩阵的对角项，其他所有需要的项都可以从A中推断出来。E的对角项可以有如下公式计算：

2D:
3D：

这些A的相应值需要根据实际情况赋予，比如对于非流体，需要用0替代。

改进的不完全乔里斯基(MIC)

我们可以在没有额外花销的情况下，获得更好的提升，MIC可以让我们需要迭代的次数方根下降。改进的主要地方是不直接丢弃我们不想要的非零值，而是考虑到他们，把他们加到对角上。

MIC（0）构造了一个下三角矩阵L，和A的下三角部分非零项一致。
1）A的非对角非零项和LLT对应位置项相等
2）A的每一行的和和LLT每一行的和相等

改进后的不完全乔里斯基的E的对角项的计算有一点不同，下面给出不同维度的公式：

因此改进后的预处理器的计算流程：

需要注意的是，这里没有直接使用改进后的不完全乔里斯基法计算预处理器，而是和没有改进的版本做了一个加权平均，只不过改进后的版本占了0.97的高权重。而且公式中为了避免除法(前面的计算E对角项公式中，计算时需要计算A/E)，我们存储了E的倒数，并且加上一个极小值，避免除0.

计算完预处理器E的值后，利用MICCG(0)算法计算压力方程的算法流程为：

Project总结

上面关于压力方程的章节的内容整合一下，可以构成压力方程的整个流程：

计算压力方程右侧的散度d，在固体边界处特殊处理
设置系数矩阵A的项，比如存储Adiag，APlusi等等
构建改进的不完全乔里斯基预处理器
用MICCG（0）算法解方程Ap=d，解得压力p
最后利用计算出来的压力根据压力梯度更新速度u。

MICCG(0)算法就是利用0级改进的不完全乔里斯基作为预处理器的预处理器共轭梯度算法

精确曲面边界

最后完善一下前面留下的大概最后一个问题，就是前面的边界都是假设和网格对齐的，而且都是完全填充的，不存在网格单元的一部分是水一部分是空气的情况，现在就来进行扩展，扩展到任意边界。

自由面的精确边界已经被很好地解决了，利用一种“幽灵”压力的思想，在液体和气体之间进行线性插值计算压力P的时候，我们取真正的液体边界上的压力为0，而不是取气体单元的中心点。

举一个一维的例子，假设边界在第i和第i+1个单元之间，比例是θ，也就是说xF = (1 − θ)xi + θxi+1处是自由面的真正位置。我们现在想要做的是k可以隐含获取 p(xF) = 0，但是不需要真的在那个位置增加一个压力值(这也会破坏网格的结构),

观察到xi+1处真正的压强是0，但是在这个自由面上不平滑，在到i之前一直非零，然后突变到0，利用这样不连续的方程不能很好地构建有限差分计算导数，所以我们考虑压强是连续变化的，从水中一直连续到空气中，并且途中在xF出为0，满足这几个条件可以构建：

线性插值检查一下在xF处确实是0.接下来就简单了，我们之前都是把相应pi的值直接设置为0，如果pi是空气网格的话。但是现在我们需要使用这个新公式来代替0：

上式是在压力更新公式中的改进，而在压力方程中，我们只需要把加到对角项，而右侧的散度部分不需要任何改变。

poi～

石星