现假设我们有一流体格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.
流体力学,此文作为科普很不错~