关于MAC网格不可压流体求解中投影矩阵的笔记

现假设我们有一流体格i，和与之相邻的格子j。我们研究投影矩阵A中，代表(i,j)格子对的元素（下简记为A[i,j]）和i格子所对应的对角元（下简记为A[i,i]）的数值情况。注意，i和j都是格子，本身坐标是一个向量，不要与矩阵坐标相混淆。

现列举j的几种情况。

情况一：j亦是流体格。此时，j格压强也是未知数，A[i,j]=-1，而A[i,i]中有一项1，如下图。（或在porous flow中，面积分数为a，那就是-a和a）

情况二：j亦是流体格，但i和j中间有一Neumann边界条件（虽然这很奇怪，我怀疑它会不会出现）。此时格子i的方程中，i-j相关联的项直接被取消，因此A[i,j]和A[i,i]的相应项都是0，如下图。

情况三：j是固体格。此时，i-j之间必有一Neumann边界条件，否则边界条件设置错误。此时，应强行将界面速度设为0，而且在矫正过程中，给j格施加一人工压强，值和i格的压强相等。从数学上，这样就使得矫正步不会改变界面处的速度。因此，i格对应的方程（也就是投影矩阵A的一行）中，i和j的压强项相互抵消，因此A[i,j]和A[i,i]中的相应项均为0，如下图。

情况四：j在计算域之外，而i和j之间没有边界条件。此时，意味着i-j之间是计算域的自由边界，一般认为j是空气，压强为常数（为方便起见，一般认为是0）.因此，i格对应的方程中，有i的压强项而无j的压强项。这意味着A[i,j]=0，A[i,i]有一项1，如下图。

情况五：j是空气格。此时，求解的是一个带自由表面的不可压流体问题。如果没有LevelSet等辅助信息追踪表面，则可简单地归结为情况三，也就是A[i,j]=0，A[i,i]有一项1.而如果有LevelSet，则可根据i，j两格的phi函数，求出来一个修正系数theta，详见：http://43.135.142.98/2021/07/11/free-surface-theta/

所以仍然A[i,j]=0，A[i,i]的那一项是那个修正系数theta，如下图。

总的规律是，当且仅当i,j都是纯流体格且中间无边界条件时，A[i,j]=-1（在porous flow中是-a），其余情况下都是0.

至于A[i,i]的这一项，若有Neumann边界条件则为0。如果是水-气界面，则为theta，否则就是1 （在porous flow中是a） .

此外还需注意一个事实：如果i和j都是流体格，则i的方程中关于j的一项，和j的方程中关于i的一项是对称的，都是1（或a）.因此我们可以在MAC网格的面节点上，把所有这样的A[i,j]存下来。这里注意，有的面可能两边都不是流体格，此时这个A[i,j]的值无意义，不出现在最后的泊松方程组中。在编程时可以随便赋个值，比如说0.