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

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

现列举j的几种情况。

情况一:j亦是流体格。此时,j格压强也是未知数,A[i,j]=-1(或在porous flow中,可能是相应的面积分数),而A[i,i]中有一项1,如下图。

情况二: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,详见:https://wmdcstdio.com/2021/07/11/free-surface-theta/

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

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

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

此外还需注意一个事实:如果i和j都是流体格,则i的方程中关于j的一项,和j的方程中关于i的一项是对称的,都是1.因此我们可以在MAC网格的面节点上,把所有这样的A[i,j]存下来。

《关于MAC网格不可压流体求解中投影矩阵的笔记》有一个想法

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注