LBGK与MRT-BGK模型计算结果对比

LBGK是单松弛模型，即f1-f9，每个离散速度方向的平衡态分布函数所采用的弛豫时间是相同的，这不是很合理，所以需要采用多松弛的模型，这个模型也很简单，就是相当于采用不同的弛豫时间，但是实际上实现起来并不是简单的给九个不同的弛豫时间，而是涉及到一个相空间转换的过程，但是归结到最后实际上是采用一个$9\times 9$的矩阵来乘以f1-f9组成的向量。

简单写一下公式
1、对于LBGK
$${\mathbf{f}}^{\ast }(\mathbf{r},t)=\mathbf{f}(\mathbf{r},t)-\frac{1}{\tau }(\mathbf{f}(\mathbf{r},t)-{\mathbf{f}}^{eq}(\mathbf{r},t))$$

2、对于MRT-BGK
$${\mathbf{f}}^{\ast }(\mathbf{r},t)=\mathbf{f}(\mathbf{r},t)-A(\mathbf{f}(\mathbf{r},t)-{\mathbf{f}}^{eq}(\mathbf{r},t))$$
A是一个$9\times 9$的矩阵
$$A=M^{-1}SM$$
M是正交矩阵，对于给定的DnQ模型，其是确定的。
S是对角松弛矩阵，其有些参数于体积粘度和动力粘度有关，有些是可调参数。

注意：上面两个公式描述的碰撞过程，LBGK与MRT-BGK仅在碰撞过程有差异。
这里其实有个疑问，我在参考资料里看到有些人说在速度空间中碰撞过程是难以实现的，在矩空间中容易实现。这句话我比较存疑，无论从公式还是代码实现来看，似乎没什么区别，甚至转换到矩空间还要进行多余的运算，所以我都是直接在速度空间进行碰撞。下面给出转换到矩空间碰撞的公式，有了解的小伙伴可以一起交流学习：
1、从速度空间转换到矩空间
$$\mathbf{m}=M\mathbf{f}$$
2、在矩空间中进行碰撞
$${\mathbf{m}}^{\ast }=S(\mathbf{m}-{\mathbf{m}}_{eq})$$
3、从矩空间转换为速度空间
$${\mathbf{f}}^{\ast }=M^{-1}{\mathbf{m}}^{\ast }$$

其实还有另一种写法，差别不大，但是这里也给出
1、从速度空间转换到矩空间
$$\mathbf{m}=M\mathbf{f}$$
2、对上述第一种方法左乘$M^{-1}$
$${\mathbf{f}}^{\ast }=M^{-1}S(\mathbf{m}-{\mathbf{m}}_{eq})$$

我个人感觉是没有区别的。
这里还补充一下${\mathbf{m}}_{eq}$的求法，在何雅玲的书中有提到，MRT方法的第二个好处是平衡态分布函数可以直接根据目标方程选取，而不是直接取$M{\mathbf{f}}_{eq}$。但是我在编程实现的过程中还是简单直接上速度空间的平衡态分布函数转换到矩空间这种方法。

OK，下面说一下MRT比LBGK计算结果的优点。
老师的PPT里有一页进行了说明

根据流线这个算例应该是Re=2000时，在方腔的左上角的涡处。

我对Re=1000的结果进行了计算，下面给出的我的结果
这个是LBGK

下面是MRT

（这里的涡量似乎大小有差异。。。）目前还没找到这个bug在哪。这两个算例的密度取值不同，但是应该不影响。。。
刚刚对LBGK的结果用tecplot计算了一下涡量，发现与程序输出的涡量确实存在差别，两个结果都用tecplot计算涡量时，结果很相近。
LBGK:wmin:-185-80wmax
MRT-LBM:wmin-162-66wmax
这里说明一下，LBGK的程序是我参考老师发的，里面有计算涡量的输出，MRT是我自己编的，没有涡量输出。
现在看来，问题应该出在tecplot计算的涡量和LBGK程序计算的涡量之间。

anyway,这不是重点。

1、LBGK

2、MRT-BGK

以上给出左上角区域放大图，可以看出MRT的涡量等值线确实光滑了许多。