摘要:为提高多松弛(MRT)格子Boltzmann模型的计算效率, 运用反演法提出了一个求解三维不可压缩流的12速MRT格子Boltzmann模型(iD3Q12 MRT模型). 这个模型比通常使用的D3Q13 MRT模型具有更高的计算效率. 在数值模拟部分我们把iD3Q12 MRT模型与可压缩性较小的一个13速多松弛模型(He-Luo D3Q13 MRT模型)在精确性和稳定性方面作比较. 通过模拟不同的流动, 包括压力驱动的稳态泊肃叶流、周期变化的压力驱动的非稳态脉动流、顶盖驱动的方腔流, 可以发现iD3Q12 MRT模型模拟以上三种流动时得到的数值解与解析解或与已有的结果符合很好, 这说明我们提出的iD3Q12 MRT模型是准确的. 在模拟稳态的泊肃叶流时, 两个模型计算的速度场的全局相对误差完全相同, 且两个模型都具有二阶的空间精度. 在模拟非稳态脉动流时, 大多情况下是12速模型的计算误差更小, 但在脉动流的最大压降增大时, iD3Q12 MRT模型先发散, 这说明He-Luo D3Q13 MRT模型具有更好的稳定性. 在模拟不同雷诺数下的顶盖驱动的方腔流时, He-Luo D3Q13 MRT模型也比iD3Q12 MRT模型更稳定.
关键词: 三维12速/
多松弛/
格子Boltzmann方法/
不可压缩流

English Abstract

全文HTML

--> --> -->

格子Boltzmann方法(LBM)起源于20世纪70年代提出的格子气自动机方法(LGA), 它克服了LGA方法的一些缺陷, 例如消除了统计噪声且其对应的宏观方程满足伽利略不变性^[1–3]. 并且LBM自身具有良好的计算局部性、程序的简洁性和拓展性等优点. 另一方面从理论上, LBM方法可以从连续Boltzmann方程得到^[4–6]. 因此, LBM作为计算流体动力学(CFD)中一种有效的介观数值模拟方法受到了广泛的关注. 它被应用在模拟一些复杂的流体, 例如: 多相流^[7–9]、悬浮液^[10]、磁流体^[11,12]、多孔介质中的流体^[13], 还被用于求解一些偏微分方程, 例如Burges方程^[14]、对流扩散方程^[15,16]、泊松方程^[17]、分数阶扩散方程^[18].
在研究人员的努力和发展下, LBM出现了lattice Bhatnagar-Gross-Krook (LBGK)模型或称为单松弛(SRT)模型^[19–23]、熵模型^[24]、双松弛(TRT)模型^[25,26]、多松弛(MRT)模型^[27–34]等. 其中LBGK模型因其形式的简洁性在各种复杂的流体传输问题研究中受到了广泛的应用, 其中最具有代表性的是Qian等^[20]提出的DdQq模型. 但由于LBGK模型以单松弛时间近似为基础, 其在稳定性和精确性方面存在不足. 几乎同时, 法国****d'Humières^[27]提出了广义的格子Boltzmann模型即多松弛格子Boltzmann模型(MRT). MRT模型和LBGK模型的不同之处主要体现在它们的碰撞项, MRT模型具有数量最多的自由度, 即可调松弛因子的数量最多, 而不同的松弛因子可以最大限度地优化模型性质, 如稳定性和精确性^[30,35], 因此MRT模型受到越来越多的重视. 常见的MRT模型有二维的D2Q9 MRT模型^[27], 三维的D3Q15 MRT模型^[28]、D3Q19 MRT模型^[28], 至今为止三维空间中速度方向最少的MRT模型是由法国****d'Humières等^[29]提出的D3Q13 MRT模型. 这些模型都可以通过Chapman-Enskog (C-E)展开恢复到Navier-Stokes方程 (N-S方程). 速度集合多的模型具有更好的稳定性和各向同性^[36,37], 但同时在计算时间和存储空间的消耗上会增加. 为此有****提出了拥有更少速度方向的MRT模型, 如D2Q8 MRT模型^[32]、D3Q14 MRT模型和D3Q18 MRT模型^[33]. 这几种MRT模型的构造原理是基于Guo等^[22]提出的不可压DdQq LBGK模型的宏观量计算与分布函数$ f_{0} $无关, 且在DdQq MRT模型^[28]恢复宏观方程的过程中没有用到能量的平方项. 因此将DdQq MRT模型的离散速度集合舍弃0速度方向、模型的矩中能量平方项丢掉, 并且将它的变换矩阵去掉第一列和能量平方项所对应的那一行再正交化得到新的DdQ (q-1) MRT模型的变换矩阵. 最后将构造的变换矩阵乘以原LBGK模型中去掉0方向的平衡态分布函数就得到DdQ (q-1) MRT模型的矩平衡态.
本文在D3Q13 MRT模型^[29]的基础上提出了求解不可压缩N-S方程的三维12速MRT模型(iD3Q12 MRT), 和已有的D2Q8, D3Q14, D3Q18等MRT模型的构造方法相比, 最大的难点在于D3Q13 MRT模型没有对应的LBGK模型, 因此就不能通过变换矩阵乘以平衡态分布函数的方法得到矩平衡态. 但可仿照D2Q8, D3Q14, D3Q18 MRT模型的构造方法构造离散速度集合和变换矩阵, 舍弃了原13速MRT模型的矩中能量e的那一项, 并且在13速MRT模型没有对应的LBGK模型的情况下运用反演法, 即用13速MRT模型变换矩阵的逆乘它的矩平衡态得到形式上的“平衡态分布函数”, 再通过一系列构造变成含有12个元素的“平衡态分布函数”, 最后将构造出的12速MRT模型的变换矩阵乘以这个“平衡态分布函数”就得到了矩平衡态, 使得新的iD3Q12 MRT模型可以在低马赫数(Ma)条件下通过C-E展开恢复到不可压N-S方程. 需要注意的是, 本文在模拟部分与iD3Q12 MRT模型作对比的是D3Q13 MRT模型的不可压版本, 这里称之为He-Luo D3Q13 MRT模型. 这是因为它的矩平衡态是将D3Q13 MRT模型的矩平衡态表达式中与速度相乘的密度ρ变成$ \rho_0 $, 这与He-Luo LBGK模型^[21]的平衡态分布函数的构造方法相同. 这种变化忽略了平衡态分布函数中Ma三次方的同阶或高阶无穷小量, 从而减少了LB模型的可压缩效应. 因此在模拟部分, 选取可压缩误差更小的He-Luo D3Q13 MRT模型与iD3Q12 MRT模型作对比.
总之, 基于D3Q13 MRT模型, 运用反演法提出了一个可以在低Ma假设下恢复到不可压N-S方程的iD3Q12 MRT模型, 它可能是目前三维MRT模型中离散速度方向个数最少的一个模型, 因此iD3Q12 MRT模型在计算量和存储量的需求上更小. 通过一系列数值模拟, 我们将iD3Q12与He-Luo D3Q13 MRT模型作对比, 验证了我们提出的iD3Q12 MRT模型的有效性, 并考察了该模型在精确性和稳定性方面与He-Luo D3Q13 MRT模型的差异.

D3Q13 MRT模型是d'Humières等^[29]提出的. 在现有的MRT模型中, 它满足伽利略不变性和各向同性, 并且是能够通过C-E展开恢复到Navier-Stokes方程的具有最少的离散速度的一个MRT模型. 在假设空间步长$ {\text{δ}{x}} = 1 $和时间步长$ {\text{δ}{t}} = 1 $, 即粒子速度$ c = {{\rm{\text{δ}}}{x}}/{{\rm{\text{δ}}}{t}} = 1 $的情形下, D3Q13 MRT模型选取的离散速度如下:

$\{{{{ c}}_0},{{{ c}}_1},\cdots,{{{ c}}_{12}}\} = \left( {\begin{array}{*{20}{c}}0&1&1&1&1&0&0&{ - 1}&{ - 1}&{ - 1}&{ - 1}&0&0\\0&1&{ - 1}&0&0&1&1&{ - 1}&1&0&0&{ - 1}&{ - 1}\\0&0&0&1&{ - 1}&1&{ - 1}&0&0&{ - 1}&1&{ - 1}&1\end{array}} \right),$

模型的演化方程为

$ \begin{split} & {f_i }({{x}} + {{{c}}_i }{\text{δ}t},t + {\text{δ}t}) - {f_i }({{x}},t) \\ =& - \sum\limits_{j = 0}^{12} {{\varLambda _{i j}}({f_j}({{x}},t) \!-\! f_j^{({\rm{eq}})}({{x}},t))}, \; i = 0\text{---}12, \end{split} $

$ f_i({{{x}}}, t) $是沿速度 $ {{{c}}}_i $移动的粒子的分布函数, $ f_i^{({\rm{eq}})}({{{x}}}, t) $是平衡态分布函数, $ \varLambda_{i j} $是$ 13\times13 $阶碰撞矩阵Λ的元素. 由速度方向$ {c_{i}} $生成了两两正交的向量$ {e_{k}}, \; k\in{0, 1, \cdots, 12} $, 再由这13个正交向量$ { e}_{k} $定义D3Q13 MRT模型的变换矩阵$ { T} $,

$ {{ T}} = \left( {\begin{array}{*{18}c} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 1 & 1 & 1 & 1 & 0 & 0 & -1 &-1 &-1 &-1 & 0 & 0 \\ 0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & 1 & 0 & 0 &-1 & -1 \\ 0 & 0 & 0 & 1 &-1 & 1 &-1 & 0 & 0 & -1 & 1 &-1 & 1 \\ -12 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 1 & 1 & 1 & 1 & -2 & -2 & 1 & 1 & 1 & 1 &-2 & -2 \\ 0 & 1 & 1 & -1 & -1 & 0 & 0 & 1 & 1 & -1 & -1 & 0 & 0 \\ 0 & 1 &-1 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 \\ 0 & 0 & 0 & 1 &-1 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & 1 & 1 &-1 &-1 & 0 & 0 & -1 & -1 & 1 & 1 & 0 & 0\\ 0 & -1 & 1 & 0 & 0 & 1 & 1 & 1 & -1 & 0 & 0 & -1 & -1 \\ 0 & 0 & 0 & 1 &-1 &-1 & 1 & 0 & 0 & -1 & 1 & 1 & -1 \end{array} } \right). $

作用于分布函数f得到13个矩$ m_{k} = {{{e}}_{k}}\cdot{{{f}}} =$$ \displaystyle\sum\limits_{i = 0}^{12}f_{i}e_{ki}, \; k\in{0, 1, \cdots, 12},$

$\begin{split}{ m}( x,t) =\,& (\rho,j_x,j_y,j_z,e,3S_{xx},S_{\omega \omega},S_{xy},\\ &S_{yz},S_{xz},h_{x},h_{y},h_{z}),\end{split}$

选取的矩平衡态$ {{m}}^{({\rm{eq}})} $为

$ \rho^{({\rm{eq}})} = \rho,\tag{4a}$

$ \left\{ \begin{aligned} & j_x^{({\rm{eq}})} = \rho{u_x}, \\ & j_y^{({\rm{eq}})} = \rho{u_y}, \\ & j_z^{({\rm{eq}})} = \rho{u_z}, \end{aligned} \right. \tag{4b}$

$ e^{({\rm{eq}})} = \frac{3}{2}(13c^2_s-8)\rho+\frac{13}{2}\rho{{{u}}^2}, \tag{4c}$

$ \left\{ \begin{aligned} & 3S_{xx}^{({\rm{eq}})} = \rho(3u_x^2-u^2), \\ & S_{\omega \omega}^{({\rm{eq}})} = \rho(u_{y}^2-u_{z}^2), \end{aligned} \right. \tag{4d}$

$ \left\{ \begin{aligned} & S_{xy}^{({\rm{eq}})} = \rho(u_x u_y), \\ & S_{yz}^{({\rm{eq}})} = \rho(u_y u_z), \\ & S_{xz}^{({\rm{eq}})} = \rho(u_x u_z), \end{aligned} \right. \tag{4e} $

$ \left\{ \begin{aligned} & h_x^{({\rm{eq}})} = 0, \\ & h_y^{({\rm{eq}})} = 0, \\ & h_z^{({\rm{eq}})} = 0, \end{aligned} \right. \tag{4f} $

通过C-E展开, 并在低Ma假设下, 该模型能够恢复到可压缩的Navier-Stokes方程:

$ {{\partial{\rho}} \over {\partial t}}+\nabla \cdot{(\rho{{u}})} = 0, \tag{5a}$

$ \begin{split} & {{\partial {(\rho{{u}}})} \over {\partial t}} +\nabla \cdot {(\rho{{u}}{{u}})} \\ =\, & - c_s^2{\nabla\rho} + {\rho}{\nu} {\Delta}{({{u}})} +{\rho} \left( {\frac{\nu }{3} + \zeta } \right)\nabla(\nabla\cdot{{{u}}}), \end{split}\tag{5b}$

其中

$\begin{split} & \nu = \frac{1}{4}\left( {\frac{1}{\lambda_{\nu}}-\frac{1}{2}} \right) = \frac{1}{2}\left( {\frac{1}{\lambda'_{\nu}}-\frac{1}{2}} \right),\\ & \zeta = \left( {\frac{2}{3}-c_{\rm s}^{2}} \right)\left( {\frac{1}{\lambda_{e}}-\frac{1}{2}} \right),\end{split} $

分别表示剪切黏度和体黏度, $ \lambda_{\nu} $和$ \lambda_{\nu}^{\prime} $是与剪切黏度ν相关的碰撞因子.
He-Luo D3Q13 MRT模型的矩平衡态选取如下:

$ \rho^{({\rm{eq}})} = \rho,\tag{7a}$

$ \left\{ \begin{aligned} & j_x^{({\rm{eq}})} = {\rho_{0}}u_x, \\ & j_y^{({\rm{eq}})} = {\rho_{0}}u_y, \\ & j_z^{({\rm{eq}})} = {\rho_{0}}u_z, \end{aligned} \right. \tag{7b}$

$ e^{({\rm{eq}})} = \frac{3}{2}(13c^2_s-8)\rho+\frac{13}{2}{{\rho}_0}{{{u}}^2},\tag{7c}$

$ \left\{ \begin{aligned} & 3S_{xx}^{({\rm{eq}})} = {\rho_{0}}(3u_x^2-u^2), \\ & S_{\omega \omega}^{({\rm{eq}})} = {\rho_{0}}(u_{y}^2-u_{z}^2), \end{aligned}\right.\tag{7d}$

$ \left\{\begin{aligned} & S_{xy}^{({\rm{eq}})} = {\rho_{0}}u_x u_y, \\ & S_{yz}^{({\rm{eq}})} = {\rho_{0}}u_y u_z, \\ & S_{xz}^{({\rm{eq}})} = {\rho_{0}}u_x u_z, \end{aligned}\right. \tag{7e} $

$ \left\{ \begin{aligned} & h_x^{({\rm{eq}})} = 0, \\ & h_y^{({\rm{eq}})} = 0, \\ & h_z^{({\rm{eq}})} = 0, \end{aligned} \right.\tag{7f}$

通过C-E展开可以将He-Luo D3Q13 MRT模型恢复到如下形式的Navier-Stokes方程:

$ \frac{{\partial}{P}}{{\partial}{t}}+{\nabla}\cdot{{{u}}} = 0,\tag{8a}$

$ {\frac{\partial{{{u}}}}{\partial{t}}}+ {\nabla}\cdot{({{{u}}}{{{u}}})} = - {\nabla {P}} + \nu {\Delta}{{{u}}}+ \left( {\frac{\nu}{3}+{\zeta}} \right) {\nabla}({\nabla}\cdot{{{u}}}),\tag{8b}$

其中

$\begin{split} &\nu = \frac{1}{4}\left( {\frac{1}{\lambda_{\nu}}-\frac{1}{2}} \right) = \frac{1}{2} \left( {\frac{1}{\lambda'_{\nu}}-\frac{1}{2}} \right),\\ & \zeta = \left( {\frac{2}{3}-c_{\rm s}^{2}} \right)\left( {\frac{1}{\lambda_{e}}-\frac{1}{2}} \right),\\ & P = p/{\rho_{0}} = c_s^2{\rho}/{{\rho}_{0}}. \end{split} $

在低Ma假设和$ {T}\gg{{L}/{c_s}} $(T和L分别表示特征时间和特征长度)的条件下, He-Luo D3Q13 MRT模型可以进一步恢复到不可压N-S方程:

$ {\nabla}\cdot{{{u}}} = 0,\tag{10a}$

$ \frac{\partial{{{u}}}}{\partial{t}}+{{{u}}}\cdot{{\nabla}{{{u}}}} = -{\nabla}{P}+{\nu}{\Delta}{{{u}}}.\tag{10b}$

在He-Luo D3Q13 MRT模型中, 宏观量的计算格式如下:

$ \rho = \sum\limits_{i = 0}^{12}f_{i},\quad \rho_{0}{{{u}}} = \sum\limits_{i = 0}^{12}{{{{c}}}_{i}}f_{i}, $

在后面的模拟中, 假设$ \rho_{0} = 1 $.

将D3Q13 MRT模型变换矩阵的逆乘它的矩平衡态得到形式上的“平衡态分布函数”, 如下所示:

$ \left(\!\!\begin{array}{*{20}c} {-(\rho*(3c_{\rm s}^2+u^2-2))/2}\\ {(\rho*(c_{\rm s}^2+u_{x}^2+2u_{x}u_{y}+u_{x}+u_{y}^2+u_{y}-u_{z}^2))/8}\\ {(\rho*(c_{\rm s}^2+u_{x}^2-2u_{x}u_{y}+u_{x}+u_{y}^2-u_{y}-u_{z}^2))/8}\\ {(\rho*(c_{\rm s}^2+u_{x}^2+2u_{x}u_{z}+u_{x}-u_{y}^2+u_{z}^2+u_{z}))/8}\\ {(\rho*(c_{\rm s}^2+u_{x}^2-2u_{x}u_{z}+u_{x}-u_{y}^2+u_{z}^2-u_{z}))/8}\\ {(\rho*(c_{\rm s}^2+u^2-2u_{x}^2+u_{y}+u_{z}+2u_{y}u_{z}))/8}\\ {(\rho*(c_{\rm s}^2+u^2-2u_{x}^2+u_{y}-u_{z}-2u_{y}u_{z}))/8}\\ {(\rho*(c_{\rm s}^2+u_{x}^2+2u_{x}u_{y}-u_{x}+u_{y}^2-u_{y}-u_{z}^2))/8}\\ {(\rho*(c_{\rm s}^2+u_{x}^2-2u_{x}u_{y}-u_{x}+u_{y}^2+u_{y}-u_{z}^2))/8}\\ {(\rho*(c_{\rm s}^2+u_{x}^2+2u_{x}u_{z}-u_{x}-u_{y}^2+u_{z}^2-u_{z}))/8}\\ {(\rho*(c_{\rm s}^2+u_{x}^2-2u_{x}u_{z}-u_{x}-u_{y}^2+u_{z}^2+u_{z}))/8}\\ {-(\rho*(-c_{\rm s}^2-u^2+2u_{x}^2+u_{y}+u_{z}-2u_{y}u_{z}))/8}\\ {-(\rho*(-c_{\rm s}^2-u^2+2u_{x}^2+u_{y}-u_{z}+2u_{y}u_{z}))/8} \end{array}\!\!\right) . $

由于He-Luo LBGK模型^[21]的平衡态分布函数的构造是基于Qian等^[20]的LBGK模型, 并将其平衡态分布函数中速度u前的密度ρ替换为$ \rho_0 $. 因此, 将(12)式中速度u前的ρ替换成$ \rho_0 $, 得到He-Luo LBGK模型式的“平衡态分布函数”, 如下所示:

$ \left(\!\!\!\begin{array}{*{20}c} {-(\rho*(3c_{\rm s}^2-2)+\rho_{0}u^2))/2}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u_{x}^2+2u_{x}u_{y}+u_{x}+u_{y}^2+u_{y}-u_{z}^2))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u_{x}^2-2u_{x}u_{y}+u_{x}+u_{y}^2-u_{y}-u_{z}^2))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u_{x}^2+2u_{x}u_{z}+u_{x}-u_{y}^2+u_{z}^2+u_{z}))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u_{x}^2-2u_{x}u_{z}+u_{x}-u_{y}^2+u_{z}^2-u_{z}))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u^2-2u_{x}^2+u_{y}+u_{z}+2u_{y}u_{z}))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u^2-2u_{x}^2+u_{y}-u_{z}-2u_{y}u_{z}))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u_{x}^2+2u_{x}u_{y}-u_{x}+u_{y}^2-u_{y}-u_{z}^2))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u_{x}^2-2u_{x}u_{y}-u_{x}+u_{y}^2+u_{y}-u_{z}^2))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u_{x}^2+2u_{x}u_{z}-u_{x}-u_{y}^2+u_{z}^2-u_{z}))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u_{x}^2-2u_{x}u_{z}-u_{x}-u_{y}^2+u_{z}^2+u_{z}))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u^2-2u_{x}^2-u_{y}-u_{z}+2u_{y}u_{z}))/8}\\ {(\rho*c_{\rm s}^2+\rho_{0}*(u^2-2u_{x}^2-u_{y}+u_{z}-2u_{y}u_{z}))/8} \end{array}\!\!\!\right), $

由于不可压LBGK模型^[22]的平衡态分布函数可以看作是在He-Luo LBGK模型的基础上, 非0方向的平衡态分布函数中的ρ换成$ p/{c_s^2} $, 与速度相乘的$ \rho_0 $令成1, 0方向的“平衡态分布函数”定义为$ \rho_0 $减去其他方向的平衡态分布函数之和. 按照这种方法, 变化(13)式, 则得到不可压LBGK模型式的“平衡态分布函数”, 如下所示:

$ \left(\!\!\!\begin{array}{*{20}c} {-(3/2)p+\rho_{0}-u^2/2}\\ {(p+(u_{x}^2+2u_{x}u_{y}+u_{x}+u_{y}^2+u_{y}-u_{z}^2))/8}\\ {(p+(u_{x}^2-2u_{x}u_{y}+u_{x}+u_{y}^2-u_{y}-u_{z}^2))/8}\\ {(p+(u_{x}^2+2u_{x}u_{z}+u_{x}-u_{y}^2+u_{z}^2+u_{z}))/8}\\ {(p+(u_{x}^2-2u_{x}u_{z}+u_{x}-u_{y}^2+u_{z}^2-u_{z}))/8}\\ {(p+(u^2-2u_{x}^2+u_{y}+u_{z}+2u_{y}u_{z}))/8}\\ {(p+(u^2-2u_{x}^2+u_{y}-u_{z}-2u_{y}u_{z}))/8}\\ {(p+(u_{x}^2+2u_{x}u_{y}-u_{x}+u_{y}^2-u_{y}-u_{z}^2))/8}\\ {(p+(u_{x}^2-2u_{x}u_{y}-u_{x}+u_{y}^2+u_{y}-u_{z}^2))/8}\\ {(p+(u_{x}^2+2u_{x}u_{z}-u_{x}-u_{y}^2+u_{z}^2-u_{z}))/8}\\ {(p+(u_{x}^2-2u_{x}u_{z}-u_{x}-u_{y}^2+u_{z}^2+u_{z}))/8}\\ {(p+(u^2-2u_{x}^2-u_{y}-u_{z}+2u_{y}u_{z}))/8}\\ {(p+(u^2-2u_{x}^2-u_{y}+u_{z}-2u_{y}u_{z}))/8} \end{array}\!\!\!\right), $

去掉(14)式中的$ f_{0} $, 得到下面1—12方向上的“平衡态分布函数”

$ \left(\!\!\!\begin{array}{*{20}c} {(p+(u_{x}^2+2u_{x}u_{y}+u_{x}+u_{y}^2+u_{y}-u_{z}^2))/8}\\ {(p+(u_{x}^2-2u_{x}u_{y}+u_{x}+u_{y}^2-u_{y}-u_{z}^2))/8}\\ {(p+(u_{x}^2+2u_{x}u_{z}+u_{x}-u_{y}^2+u_{z}^2+u_{z}))/8}\\ {(p+(u_{x}^2-2u_{x}u_{z}+u_{x}-u_{y}^2+u_{z}^2-u_{z}))/8}\\ {(p+(u^2-2u_{x}^2+u_{y}+u_{z}+2u_{y}u_{z}))/8}\\ {(p+(u^2-2u_{x}^2+u_{y}-u_{z}-2u_{y}u_{z}))/8}\\ {(p+(u_{x}^2+2u_{x}u_{y}-u_{x}+u_{y}^2-u_{y}-u_{z}^2))/8}\\ {(p+(u_{x}^2-2u_{x}u_{y}-u_{x}+u_{y}^2+u_{y}-u_{z}^2))/8}\\ {(p+(u_{x}^2+2u_{x}u_{z}-u_{x}-u_{y}^2+u_{z}^2-u_{z}))/8}\\ {(p+(u_{x}^2-2u_{x}u_{z}-u_{x}-u_{y}^2+u_{z}^2+u_{z}))/8}\\ {(p+(u^2-2u_{x}^2-u_{y}-u_{z}+2u_{y}u_{z}))/8}\\ {(p+(u^2-2u_{x}^2-u_{y}+u_{z}-2u_{y}u_{z}))/8} \end{array}\!\!\!\right). $

根据d'Humières等^[29]提出的最原始的D3Q13 MRT模型, 我们构造出了新的三维12速多松弛格子Boltzmann模型(iD3Q12 MRT模型). 首先将13速模型的0速度方向舍弃得到了12速模型的离散速度方向的集合:

$\{{{{ c}}_1},{{{ c}}_2},\ldots,{{{ c}}_{12}}\} = \left( {\begin{array}{*{20}{c}} 1 & 1 & 1 & 1 & 0 & 0 & -1 &-1 &-1 &-1 & 0 & 0 \\ 1 & -1 & 0 & 0 & 1 & 1 & -1 & 1 & 0 & 0 &-1 &-1 \\ 0 & 0 & 1 &-1 & 1 &-1 & 0 & 0 &-1 & 1 &-1 & 1 \end{array}} \right)c,$

不失一般性, 这里假设粒子速度$ c = 1 $. iD3Q12 MRT模型的演化方程为

$\begin{split}& {f_i }({{x}} + {{{c}}_i }{\text{δ} t},t + {\text{δ} t}) - {f_i }({{x}},t) \\ =\, & - \sum\limits_{j = 1}^{12} {{\varLambda _{i j}}({f_j}({{x}},t) - f_j^{({\rm{eq}})}({{x}},t))} , \\ & i = 1\text{---} 12,\end{split}$

$ f_{i}({{x}}, t) $是沿速度 $ {{c}}_i $移动的粒子的分布函数, $ f_i^{({\rm{eq}})}({{x}}, t) $是平衡态分布函数, $ \varLambda_{i j} $是 12 × 12 阶碰撞矩阵Λ的元素. 同时, 演化方程也可以写作向量形式:

$\begin{split} & |f({{x}} + {{{c}}_i }{\text{δ}t},t + {\text{δ}t})\rangle - |f({{x}},t)\rangle \\ =\, & - \varLambda (|f({{x}},t)\rangle - |{f^{({\rm{eq}})}}({{x}},t)\rangle), \end{split} $

其中$ |f({{x}}, t)\rangle = (f_1({{x}}, t), f_2({{x}}, t), \cdots, f_{12}({{x}}, t))' $是列向量, 符号′代表转置算子. 构造一个MRT模型, 变换矩阵和矩平衡态的构造至关重要. 按照构造D2Q8^[32], D3Q14和D3Q18 MRT模型^[33]的变换矩阵的方法构造了iD3Q12 MRT模型的变换矩阵, 在构造的过程中, 舍弃了13速MRT模型的矩中能量那一项, 得到的iD3Q12 MRT模型的变换矩阵如下:

$ { T}=\left( {\begin{array}{*{18}c} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 0 & 0 & -1 &-1 &-1 &-1 & 0 & 0 \\ 1 & -1 & 0 & 0 & 1 & 1 & -1 & 1 & 0 & 0 &-1 & -1 \\ 0 & 0 & 1 &-1 & 1 &-1 & 0 & 0 &-1 & 1 &-1 & 1 \\ 1 & 1 & 1 & 1 & -2 & -2 & 1 & 1 & 1 & 1 &-2 & -2 \\ 1 & 1 & -1& -1 & 0 & 0 & 1 & 1 & - 1& - 1 & 0 & 0 \\ 1 &-1 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 \\ 0 & 0 & 1 &-1 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 \\ 1 & 1 &-1 &-1 & 0 & 0 & -1 & -1 & 1 & 1 & 0 & 0\\ -1 & 1 & 0 & 0 & 1 & 1 & 1 & -1 & 0 & 0 & -1 & -1 \\ 0 & 0 & 1 &-1 &-1 & 1 & 0 & 0 & -1 & 1 & 1 & -1 \\ \end{array} } \right). $

通过变换矩阵T将分布函数f变换为矩阵m

$ \begin{split}{ m}( x,t) =\, & (P_{1},{u}_x,{ u}_y,{ u}_z,3S_{xx},S_{\omega \omega},\\ & S_{xy},S_{yz},S_{xz},h_{x},h_{y},h_{z})',\end{split} $

碰撞矩阵$ \varLambda = {{ T}}^{-1} {\hat{ \varLambda}} {{ T}} $,

$ {\hat{{\varLambda}}}\equiv {\rm diag}(\lambda_c, \lambda_c, \lambda_c, \lambda_c, \lambda_{\nu}, \lambda_{\nu}, \lambda'_{\nu}, \lambda'_{\nu}, \lambda'_{\nu}, \lambda_{t}, \lambda_{t}, \lambda_{t}), $

$ \lambda_{\rm c} $是与守恒矩相关的松弛因子, 由于守恒矩对应的平衡态是它本身, 因此这些松弛因子的取值不影响粒子的碰撞过程, 可以将它们取值为0. $ \lambda_{\nu} $, $ \lambda'_{\nu} $和$ \lambda_t $是与非守恒矩相匹配的松弛因子, 为了保持模型的稳定性, 取值一般在(0, 2). 我们选取的矩平衡态如下:

$ {P_{1}}^{({\rm{eq}})} = {|{{u}}|}^2/2+3p/2,\tag{21a} $

$ \left\{\begin{aligned} & j_x^{({\rm{eq}})} = u_x, \\ & j_y^{({\rm{eq}})} = u_y, \\ & j_z^{({\rm{eq}})} = u_z, \end{aligned} \right. \tag{21b}$

$ \left\{ \begin{aligned} & 3S_{xx}^{({\rm{eq}})} = 3u_x^2-u^2, \\ & S_{\omega \omega}^{({\rm{eq}})} = u_{y}^2-u_{z}^2, \end{aligned}\right. \tag{21c}$

$ \left\{\begin{aligned} & S_{xy}^{({\rm{eq}})} = u_x u_y, \\ & S_{yz}^{({\rm{eq}})} = u_y u_z, \\ & S_{xz}^{({\rm{eq}})} = u_x u_z, \end{aligned}\right . \tag{21d}$

$ \left\{\begin{split}& h_x^{({\rm{eq}})} = 0, \\ & h_y^{({\rm{eq}})} = 0, \\ & h_z^{({\rm{eq}})} = 0, \end{split} \right. \tag{21e}$

p是压力, $ {|{{u}}|}^2 $是速度$ {{u}} = (u_x, u_y, u_z) $模的平方. 演化方程在矩空间的形式如下:

$\begin{split} & |f({{{x}}} + {{{{c}}}_\alpha }{\text{δ}t},t + {\text{δ}t})\rangle \\ =\, & {{{ T}}^{ - 1}} [{ m}({{{x}}},t) - {\hat{ \varLambda}}({ m}( x,t) - {{ m}^{({\rm{eq}})}}( x,t))], \end{split} $

其中$ {\hat \varLambda} = {{ T}}{{\varLambda}}{{{ T}}^{ - 1}} $, $ {{ m}^{({\rm{eq}})}}( x, t) $是由(21)式定义的矩平衡态. 这样通过C-E展开, 可以将iD3Q12 MRT模型恢复到如下形式的N-S方程:

$ \frac{{\partial}{\left( {\dfrac{3p}{2}+\dfrac{|{{{u}}}|^2}{2}} \right)}} {{\partial}{t}}+{\nabla}\cdot{{{u}}} = 0,\tag{23a}$

$ {\frac{\partial{{{u}}}}{\partial{t}}}+ {\nabla}\cdot{({{{u}}}{{{u}}})} = -{\nabla {p}} + \nu {\Delta}{{{u}}}+{\xi}\nabla({\nabla}\cdot{{{u}}}),\tag{23b}$

其中

$ \begin{split} & \nu = \frac{1}{4}\left( { \frac{1}{{\lambda}_{\nu}}- \frac{1}{2}} \right){\text{δ}t} = \frac{1}{2}\left( {\frac{1}{ {\lambda}_{\nu}'}-\frac{1}{2}} \right){\text{δ} t},\\ & \xi = \frac{1}{12}\left( {{\frac{1}{{\lambda}_{\nu}}-\frac{1}{2}}} \right){\text{δ} t}.\end{split}$

若令$ \tau = {1}/{{\lambda}_{\nu}} $, 则

$ \nu = \dfrac{1}{4}\left( { \tau \!-\! \dfrac{1}{2}} \right){\text{δ}t}, ~ \xi \!=\! \dfrac{1}{12}\left( {\tau \!-\! \dfrac{1}{2}} \right){\text{δ}t}.$

在低Ma条件下, 方程(23)可以写成如下形式:

$ \nabla \cdot{{{u}}} = 0, \tag{25a}$

$ {{\partial {{{u}}}} \over {\partial t}} + {{{u}}}\cdot\nabla {{{u}}} = -{\nabla p} + \nu {\Delta}{{{u}}}.\tag{25b}$

在iD3Q12 MRT模型中, 宏观量的计算如下:

$ p = \frac{2}{3}\left(\sum\limits_{i = 1}^{12}f_{i}-\frac{{|{{u}}|}^2}{2}\right), \quad {{{u}}} = \sum\limits_{i = 1}^{12}{{{c}}_{i}}f_{i}. $

为了验证提出的iD3Q12 MRT模型的有效性, 我们模拟了三维泊肃叶流、脉动流与顶盖驱动的方腔流. 值得注意的是, 在模型的建立和推导过程中, 假设$ c ={{\rm{\text{δ}}} x}/{{\rm{\text{δ}}} t} = 1 $. 但在实际的计算中, c有时不取1. 这种情况下, 我们将宏观速度u、压力p进行单位化$ {{{u}}}' = {{{u}}}/c $, $ {p'} = p/{c^{2}} $, 再将单位化后的量代入我们提出的模型进行计算, 计算结束后再将所得到的$ {{{u}}}', p' $分别乘以c, $ c^{2} $, 就得到我们所要求的宏观速度和压力. 同时, 剪切黏性系数按照$ \nu = \dfrac{1}{4}\left( {\dfrac{1}{{\lambda}_{\nu}}- \dfrac{1}{2}} \right){c^{2}}{\text{δ} t} $确定. 此外, 碰撞矩阵$ {\hat{ \varLambda}} $中的松弛因子取定如下: $ {\lambda}_{c} = 1.0 $, $ {\lambda}_{\nu} = {1}/{\tau} $, $ {\lambda}_{t}=1.8. $我们也给出了He-Luo D3Q13 MRT模型计算得到的有关结果. 对于边界条件, 使用的是非平衡态外推法^[38]. 程序的计算流程如下:
1) 初始化　输入计算参数, 宏观速度u和压力p的初始值都设置为0, 并初始化分布函数f;
2)矩空间碰撞

$ \begin{split} \, &{{{f}}}^{+}({{x}},t) = {{f}}({{x}},t)-{\varLambda}[{{f}}({{x}},t) -{{f}}^{({\rm{eq}})}({{x}},t)] \\ =\,& {{{{ T}}}}^{-1}\{{{m}}({{x}},t)-\hat{{\varLambda}}[{{m}}({{x}},t) -{{m}}^{({\rm{eq}})}({{x}},t)]\}; \end{split} $

3) 速度空间迁移

$ {{f}}({{x}}+{{c}}_{i}\text{δ}{t},t+\text{δ}{t}) = {{{f}}}^{+}({{x}},t); $

4) 边界处理采用非平衡态外推法.
2

-->

5.1.三维稳态的泊肃叶流

如图1所示, 三维泊肃叶流的物理空间限制在一个长方体通道中, 其长、宽、高分别为$ 0\leqslant x \leqslant l $, $ -a\leqslant y \leqslant a $, $ -b\leqslant z \leqslant b $, $ l = 2 $, $ a = b = 0.5 $, 原点O代表流体入口平面的中心.
图 1 三维泊肃叶流示意图
Figure1. The schematic of three-dimensional Poiseuille flow.

边界条件设置为

$ { u}(x,\pm a,z,t) = { u}(x,y,\pm b,t) = 0, \tag{29a}$

$ p(0,y,z,t) = {p_{\rm in}},\;\;\; p(l,y,z,t) = {p_{\rm out}}, \tag{29b}$

$ p_{\rm in} $和$ p_{\rm out} $分别表示进出口压力, 三维泊肃叶流具有稳态的解析解^[39]

$ \begin{split} u_{x}(y,z,t) =\, &{{16{a^2}} \over {\nu {{\text{π}} ^3}}}\left( { - \frac{{{\rm{d}}p}}{{{\rm{d}}x}}} \right)\sum\limits_{i = 1,3,5,\cdots}^\infty {{{( - 1)}^{(i - 1)/2}}} \\& \times\left[ {1 - {{{\rm cosh}({\rm i}{\text{π}} z/(2a))} \over {{\rm cosh}({\rm i}{\text{π}} b/(2a))}}} \right] \\ & \times {{\cos({\rm i}{\text{π}} y/(2a))} \over {{i^3}}},\end{split}\tag{30a}$

$ u_{y}(x,y,z,t) = u_{z}(x,y,z,t) = 0,\tag{30b}$

$ p(x,t) = {p_{\rm in}} + {{{\rm d}p} \over {{\rm d}x}}x,\tag{30c}$

$ {\rm d}p/{\rm d}x = (p_{\rm out}-p_{\rm in})/l $是通道内的压力梯度, ν是流体的剪切黏度. 在模拟中, 初始状态是:

${ u}(x,y,z,0) = { 0}, \;\;\; p(0 < x < l,y,z,0) = {{p_{\rm in}+p_{\rm out}}\over {2}}, $

并且设$ p_{\rm in} = 1.1 $和$ p_{\rm out} = 1.0 $, 模拟达到稳定状态的判定准则是

$ {{\displaystyle\sum\nolimits_{i} | {u_{x}}( x_i,t +{\text{δ}}t) - {u_{x}}( x_i,t)|} \over {\displaystyle\sum\nolimits_{i} | {u_{x}}({{{x}}}_i,t)|}} \leqslant 1.0 \times {10^{ - 10}}, $

$\displaystyle\sum\nolimits_{i} $表示求和遍布每一个网格点. 模拟首先在网格数为$ 65\times33\times33 $的网格下进行, 且参数$ \nu = 0.03 $, $ \lambda_{\nu} = 1.3 $. 图2显示了在$ x = 1 $截面处z取不同值时水平速度$ u_{x} $随y变化的函数图像和在$ z = 0 $的截面处y取不同值时压力p随x变化的函数图像. 从图2可以看出iD3Q12 MRT模型在模拟稳态的泊肃叶流时所得到的数值解与已有的解析解符合得很好.
图 2 泊肃叶流数值解与解析解的对比　(a) 泊肃叶流在$x=1$截面处z取不同的值时水平速度$u_{x}$随y变化的函数图像; (b) 在截面$z=0$处y取不同的值时压力p随x变化的函数图像; 直线: 解析解; 符号: 数值解; 松弛因子$\lambda_{\nu}=1.3$
Figure2. Comparison between numerical and analytical solutions of Poiseuille flow: (a) The variation of $u_{x}$ with y for different locations of z at section $x=1$ for Poiseuille flow; (b) the variation of pressure with x for different locations of y at section $z=0$ for Poiseuille flow. Lines, analytical solutions; symbols, numerical results; the relaxation parameter ${\lambda}_{\nu}=1.3$.

我们也计算了iD3Q12 MRT模型在模拟泊肃叶流时的空间精度的阶. 为此, 计算了不同空间步长下的泊肃叶流的速度场的全局相对误差$ {\rm GRE}_u$, $ {\rm GRE}_u$的表达式为

$\begin{split} &{\rm {GRE}}_u \\=\, & {\sqrt{\displaystyle\sum\nolimits_i[{(u_{xn}\!-\!u_{xa})^2\!+\!(u_{yn}\!-\!u_{ya})^2\!+\!(u_{zn}\!-\!u_{za})^2}]}\over {\sqrt{\displaystyle\sum\nolimits_i[{u_{xa}^2+u_{ya}^2+u_{za}^2}]}}}, \end{split}$

公式中的n和a分别表示数值解和解析解, $ u_{x} $, $ u_{y} $和$ u_{z} $是流体速度u的三个分量, 同样的关于i的求和遍及所有网格点. 在不同松弛参数和不同空间步长下由iD3Q12 MRT和D3Q13 MRT模型计算的泊肃叶流的速度场的全局相对误差$ \rm {GRE}_{\rm u} $列在表1中. 值得注意的是, 在同样的参数下, 这两个模型的全局相对误差是完全一样的.

${\rm GRE}_u$	Lattice spacing $\text{δ} x$				Model
	1/8	1/16	1/32	1/64
${\lambda}_{\nu}=0.8,$ ${\lambda}'_{\nu}=1.143$	$3.090\times10^{-2}$	$7.700\times10^{-3}$	$1.900\times10^{-3}$	$4.623\times10^{-4}$	iD3Q12 MRT
	$3.090\times10^{-2}$	$7.700\times10^{-3}$	$1.900\times10^{-3}$	$4.623\times10^{-4}$	D3Q13 MRT
${\lambda}_{\nu}=1.0,$ ${\lambda}'_{\nu}=1.333$	$5.990\times10^{-2}$	$1.660\times10^{-2}$	$4.400\times10^{-3}$	$1.100\times10^{-3}$	iD3Q12 MRT
	$5.990\times10^{-2}$	$1.660\times10^{-2}$	$4.400\times10^{-3}$	$1.100\times10^{-3}$	D3Q13 MRT
${\lambda}_{\nu}=1.3,$ ${\lambda}'_{\nu}=1.576$	$8.720\times10^{-2}$	$2.500\times10^{-2}$	$6.700\times10^{-3}$	$1.700\times10^{-3}$	iD3Q12 MRT
	$8.720\times10^{-2}$	$2.500\times10^{-2}$	$6.700\times10^{-3}$	$1.700\times10^{-3}$	D3Q13 MRT

表1iD3Q12 MRT和D3Q13 MRT模型在不同松弛因子${\lambda}_{\nu}$和不同空间步长下计算得到的泊肃叶流的速度场的全局相对误差${\rm GRE}_u$
Table1.The ${\rm GRE}_u$ of velocity field for Poiseuille flow computed by iD3Q12 MRT and D3Q13 MRT models under different relaxation parameters and different lattice spacings.

假设$ {\rm {GRE}}_u = a(\text{δ} x)^b (a > 0) $, 然后我们得到$ \ln({\rm {GRE}}_{\rm u}) = b\ln(\text{δ} x)+\ln(a) $. 如果$ b = 2 $, 就称数值解具有二阶精度. 用表1中的数据进行最小二乘法的线性拟合, 线性拟合的图像在图3中显示. 图中三条直线分别对应着iD3Q12 MRT模型在$ \lambda_{\nu} = 0.8, 1.0 $和1.3下的拟合直线, 它们对应的斜率分别为1.94, 1.91和1.88, 都接近2. 这说明我们提出的iD3Q12 MRT模型在模拟稳态的泊肃叶流时能够达到二阶的空间精度.
图 3 不同的$\lambda_{\nu}$下, 模拟泊肃叶流得到的速度场的全局相对误差${\rm {GRE}}_u$随空间步长$\text{δ}{x}$的变化, 符号代表数值解, 连线表示拟合直线
Figure3. The variation of ${\rm {GRE}}_u$ of velocity field with the lattice spacing $\text{δ}{x}$ at different $\lambda_{\nu}$ for Poiseuille flow. Symbols represent numerical solutions, lines represent fitting line.

2 -->

5.2.三维非稳态的脉动流

模拟三维脉动流是为了验证所提出的iD3Q12 MRT模型在模拟非稳态流时的精度. 脉动流的流域与泊肃叶流的流域相同, 如图1所示. 但脉动流在管道进出口两端的压力梯度是周期性变化的, 压力梯度为

$ {{{\rm d}p}\over{{\rm d}x}} = G \cos(\omega t), $

其中G是振幅, ω是频率. 脉动流的解析解是^[40]

$\begin{split} & {u_{x}}(y,z,t) \\=\, & {\rm{Re}}\left({\rm i}\frac{G}{\omega }\left\{ 1 - 2\sum\limits_{n = 0}^\infty {\frac{{{{( - 1)}^n}}}{{{p_n}}}} \left[ \frac{{\cosh ({\gamma _n}y/b)\cos ({p_n}z/b)}}{{\cosh ({\gamma _n}a/b)}} \right.\right.\right.\\ & \left. \left. \left.+ \frac{{\cosh ({\sigma _n}z/b)\cos ({q_n}y/b)}}{{\cosh ({\sigma _n})}}\right] \right\}{\rm e}^{{\rm i}\omega t}\right),\\[-20pt]\end{split}$

其中

$ {\gamma _n} = \sqrt {p_n^2 + {\rm i}{\eta ^2}},\;\;\;\;{\sigma _n} = \sqrt {q_n^2 + {\rm i}{\eta ^2}}, \tag{36a} $

$ {p_n} = {{2n + 1} \over 2}{\text{π}},\;\;\;\; {q_n} = {{2n + 1} \over 2}{\text{π}} {b \over a}, \tag{36b}$

$ \eta = b\sqrt {\omega /\nu } $是沃门斯里(Womersley number)数. 我们将参数设置如下, 网格大小为$ 81\times41\times41 $, 压力变化的周期$ T = 100 $, 频率$ \omega = {2{\text{π}}}/{T} $, 模拟中进出口压力分别设置为$ p_{\rm in} = 0.999 $, $ p_{\rm out} = 1.0 $, 沿着管道的压力梯度的振幅$ G = {\Delta{p}}/{l} = 0.0005 $, 时间步长$ {\rm{\text{δ}}}{t} = 0.0125 $. 流动的初始状态的设置与泊肃叶流的初始状态的设置是一致的, 如(31)式所示. 计算的收敛准则是

$ {{\displaystyle\sum\nolimits_i {| {{u_x}({ x_i},t + T) - {u_x}({ x_i},t)} |} } \over {\displaystyle\sum\nolimits_i {| {{u_x}({ x_i},t + T)} |} }} \leqslant 6.0\times10^{-14}, $

关于i的求和遍及所有网格点. 图4显示了在四个不同时刻$ t = T/4, T/2, 3 T/4 $和T下, 在直线$ x = 1 $, $ z = 0 $上水平速度$ u_{x} $随y变化的函数图像, 此时$ \eta = 2.8285, \;\lambda_{\nu} = 1.522 $. 水平速度$ u_{x} $由$ U_{\max} $无量纲化, $ U_{\rm max} = 1.876\times10^{-2} $是泊肃叶流在进出口压差为$ p_{\rm out}-p_{\rm in} = -G*l $时的最大水平速度. 从图4可以得出iD3Q12 MRT模型在模拟非稳态流时得到的数值解与解析解符合较好.
图 4 在$\eta=2.8285$时脉动流在$x=1$, $z=0$处水平速度u_x随y变化的函数. 直线: 解析解; 符号: 数值解
Figure4. The variation of horizontal velocity u_x with y for pulsatile flow at the location $x=1$, $z=0$, $\eta=2.8285$. Line, analytical solutions; symbols, numerical solutions.

表2给出了在$ \eta = 2.8285 $且$ \tau = {1}/{\lambda_{\nu}} $、ν一定的条件下, 用iD3Q12 MRT和D3Q13 MRT模型模拟脉动流时得到的速度场的全局相对误差$ {\rm GRE}_u. $表3显示了由表2中的数据计算得到的相邻的两个空间步长下的iD3Q12 MRT和D3Q13 MRT模型的空间精度的阶$ \dfrac{[\ln({\rm {GRE}}_u)]_f-[\ln({\rm {GRE}}_u)]_c}{[\ln({\rm{\text{δ}}} x)]_f-[\ln({\rm{\text{δ}}} x)]_c} $, 这里f表示小的空间步长, c表示大的空间步长. 从表2中可以看出, 在四个不同时刻下由iD3Q12 MRT模型计算的速度场的全局相对误差都比D3Q13 MRT模型计算的速度场的全局相对误差要略小一些. 从表3中可以发现D3Q13 MRT模型比iD3Q12 MRT模型的阶稍微更接近于2. 根据表2中的数据, 我们绘制了用iD3Q12 MRT模型计算得到的速度场的全局相对误差$ {\rm GRE}_u $随空间步长$ \text{δ} x $变化的函数图像, 如图5所示. 在四个不同时刻t = T/4, T/2, 3T/4和T下拟合直线的斜率分别是1.98, 1.88, 1.97和1.87, 这说明iD3Q12 MRT模型在模拟非稳态脉动流时具有二阶的空间精度.

Lattice spacing	${\rm GRE}_u$				Model
	$T/4$	$T/2$	$3 T/4$	T
${\rm{\text{δ}} } x= {1}/{20}$	$1.483\times10^{-2}$	$4.214\times10^{-2}$	$1.805\times10^{-2}$	$4.028\times10^{-2}$	iD3Q12 MRT
	$1.662\times10^{-2}$	$4.733\times10^{-2}$	$2.118\times10^{-2}$	$4.299\times10^{-2}$	D3Q13 MRT
${\rm{\text{δ}} } x= {1}/{40}$	$3.803\times10^{-3}$	$1.199\times10^{-2}$	$4.651\times10^{-3}$	$1.153\times10^{-2}$	iD3Q12 MRT
	$4.172\times10^{-3}$	$1.324\times10^{-2}$	$5.398\times10^{-3}$	$1.217\times10^{-2}$	D3Q13 MRT
${\rm{\text{δ}} } x= {1}/{60}$	$1.702\times10^{-3}$	$5.569\times10^{-3}$	$2.085\times10^{-3}$	$5.369\times10^{-3}$	iD3Q12 MRT
	$1.855\times10^{-3}$	$6.116\times10^{-3}$	$2.412\times10^{-3}$	$5.648\times10^{-3}$	D3Q13 MRT
${\rm{\text{δ}} } x= {1}/{80}$	$9.605\times10^{-4}$	$3.204\times10^{-3}$	$1.177\times10^{-3}$	$3.092\times10^{-3}$	iD3Q12 MRT
	$1.043\times10^{-3}$	$3.509\times10^{-3}$	$1.360\times10^{-3}$	$3.247\times10^{-3}$	D3Q13 MRT

表2在$\eta=2.8285$时, 不同空间步长下用iD3Q12 MRT模型和D3Q13 MRT模型模拟脉动流所得的不同时刻下的速度场的全局相对误差${\rm GRE}_u$
Table2.The global relative errors of the velocity field at different times for pulsatile flow simulated by iD3Q12 MRT and D3Q13 MRT models at different lattice spacings, $\eta=2.8285$.

Adjacent spacing	Order				Model
	$T/4$	$T/2$	$3 T/4$	T
Average	1.978	1.875	1.974	1.869	iD3Q12 MRT
	1.998	1.891	1.984	1.879	D3Q13 MRT
${1}/{20} \to {1}/{40}$	1.963	1.813	1.956	1.805	iD3Q12 MRT
	1.994	1.838	1.972	1.821	D3Q13 MRT
${1}/{40}\to {1}/{60}$	1.983	1.891	1.979	1.885	iD3Q12 MRT
	1.999	1.905	1.987	1.893	D3Q13 MRT
${1}/{60} \to {1}/{80}$	1.989	1.922	1.988	1.918	iD3Q12 MRT
	2.001	1.931	1.992	1.924	D3Q13 MRT

表3相邻空间步长下的iD3Q12 MRT和D3Q13 MRT模型的空间精度的阶
Table3.The orders of the spatial accuracy of iD3Q12 MRT and D3Q13 MRT models under adjacent spacings.

图 5 同一周期四个不同时刻下变量${\rm {GRE}}_u$随空间步长$\text{δ}{x}$的变化
Figure5. The variation of ${\rm {GRE}}_u$ with the lattice spacing at four different times in a period for pulsatile flow.

当最大的压差$ \Delta{p} $增大时, LB模型计算的速度场的可压缩性误差会增大, 在这种情况下比较iD3Q12 MRT和D3Q13 MRT模型计算的速度场的全局相对误差$ \rm {GRE}_u$能凸显两个模型在精度上的差异. 表4和表5分别显示了最大的压差$ \Delta{p} $增大时iD3Q12 MRT和D3Q13 MRT模型在不同的空间步长和松弛时间τ下所计算得到的T时刻的速度场的全局相对误差. 从表4可以看出, 由iD3Q12 MRT模型计算的脉动流在T时刻的全局相对误差总是比D3Q13 MRT模型要略小一些. 这说明iD3Q12 MRT模型的计算精度比D3Q13 MRT模型的计算精度要高. 随着最大压差$ \Delta{p} $的增加iD3Q12 MRT模型的计算开始发散但D3Q13 MRT模型却没有发散, 如$ {\rm{\text{δ}}}{x} = {1}/{60} $且$ \Delta{p} = 0.08 $时. 当空间步长变为$ {\rm{\text{δ}}}{x} = {1}/{80} $时再继续增大最大压差$ \Delta{p} $同样会出现类似的情形, 如$ \Delta{p} = 0.10 $和$ \Delta{p} = $ 0.12, 这说明D3Q13 MRT模型在稳定性方面较iD3Q12 MRT模型更好. 从表5的数据可以看出不同的τ下, iD3Q12 MRT模型计算得到的速度场的全局相对误差在绝大部分情况下都比D3Q13 MRT模型的要略小一些, 但也有个别反常情形, 如$ \tau = 0.70 $且$ \Delta{p} = 0.01 $或$ \Delta{p} = 0.04 $时. 在表5中还发现和表4相似的情形, 即在$ \tau = 0.55 $时增大最大压差$ \Delta{p} $出现了iD3Q12 MRT模型先发散的情形, 如$ \Delta{p} = 0.03 $时. 从表5中还可以发现, 当τ更小时, 增加$ \Delta p $, iD3Q12 MRT和D3Q13 MRT模型容易发散. 因此, 在$ \Delta p $较大时, 适当地增大τ可以提升两个模型的稳定性.

$ \Delta p $	Lattice spacing ${\rm{\text{δ}} } x$				Model
	1/20	1/40	1/60	1/80
$0.005 $	$9.919\times10^{-2}$	$3.030\times10^{-2}$	$1.442\times10^{-2}$	$8.402\times10^{-3}$	iD3Q12 MRT
	$1.121\times10^{-1}$	$3.326\times10^{-2}$	$1.568\times10^{-2}$	$9.084\times10^{-3}$	D3Q13 MRT
$0.010$	$1.172\times10^{-1}$	$3.445\times10^{-2}$	$1.618\times10^{-2}$	$9.362\times10^{-3}$	iD3Q12 MRT
	$1.679\times10^{-1}$	$4.763\times10^{-2}$	$2.199\times10^{-2}$	$1.260\times10^{-2}$	D3Q13 MRT
$0.020$	$1.777\times10^{-1}$	$5.110\times10^{-2}$	$2.365\times10^{-2}$	$1.355\times10^{-2}$	iD3Q12 MRT
	$2.940\times10^{-1}$	$8.630\times10^{-2}$	$3.987\times10^{-2}$	$2.279\times10^{-2}$	D3Q13 MRT
$0.050$		$1.243\times10^{-1}$	$5.848\times10^{-2}$	$3.386\times10^{-2}$	iD3Q12 MRT
		$2.025\times10^{-1}$	$9.868\times10^{-2}$	$5.757\times10^{-2}$	D3Q13 MRT
$0.080$				$6.073\times10^{-2}$	iD3Q12 MRT
			$1.575\times10^{-2}$	$9.405\times10^{-2}$	D3Q13 MRT
$0.100$					iD3Q12 MRT
				$1.192\times10^{-1}$	D3Q13 MRT
$0.120$					iD3Q12 MRT
				$1.454\times10^{-1}$	D3Q13 MRT

表4在$\tau=0.5667$, $\eta=4.3416$, 最大压差$\Delta{p}$增大时不同的空间步长下由iD3Q12 MRT和D3Q13 MRT模型模拟的脉动流在时刻T下的速度场所计算的全局相对误差${\rm GRE}_u$, 空白处表示计算发散
Table4.The global relative error calculated by the velocity field at time T of pulsatile flow simulated by the iD3Q12 MRT and D3Q13 MRT models under different lattice spacings. The maximal pressure drop $ \Delta{p} $ of the channel increases, $\tau=0.5567$, $\eta=4.3416$ are fixed. The blank indicates that the computation is divergent.

$\Delta p$	τ				Model
	0.55	0.60	0.70	0.90
$0.005 $	$1.302\times10^{-1}$	$6.311\times10^{-2}$	$2.955\times10^{-2}$	$1.744\times10^{-2}$	iD3Q12 MRT
	$1.556\times10^{-1}$	$6.560\times10^{-3}$	$3.023\times10^{-2}$	$1.993\times10^{-2}$	D3Q13 MRT
$0.010$	$1.612\times10^{-1}$	$6.830\times10^{-2}$	$2.711\times10^{-2}$	$1.736\times10^{-2}$	iD3Q12 MRT
	$2.435\times10^{-1}$	$8.735\times10^{-2}$	$2.661\times10^{-2}$	$2.058\times10^{-2}$	D3Q13 MRT
$0.020$	$2.475\times10^{-1}$	$9.926\times10^{-2}$	$2.624\times10^{-2}$	$1.656\times10^{-2}$	iD3Q12 MRT
	$4.182\times10^{-1}$	$1.542\times10^{-1}$	$2.757\times10^{-2}$	$2.195\times10^{-2}$	D3Q13 MRT
$0.030$		$1.430\times10^{-1}$	$3.421\times10^{-2}$	$1.509\times10^{-2}$	iD3Q12 MRT
	$5.482\times10^{-1}$	$2.193\times10^{-1}$	$3.616\times10^{-2}$	$2.343\times10^{-2}$	D3Q13 MRT
$0.040$			$5.001\times10^{-2}$	$1.349\times10^{-2}$	iD3Q12 MRT
			$4.693\times10^{-2}$	$2.502\times10^{-2}$	D3Q13 MRT
$0.050$				$1.291\times10^{-2}$	iD3Q12 MRT
				$2.674\times10^{-2}$	D3Q13 MRT

表5在${\rm{\text{δ}}} {x}={1}/{20}$时, 最大压差$\Delta{p}$增大时不同的松弛时间τ下由iD3Q12 MRT和D3Q13 MRT模型模拟的脉动流由T时刻的速度场计算得出的全局相对误差${\rm GRE}_u$, 空白处表示计算发散
Table5.The global relative error of the velocity field at time T of the pulsatile flow simulated by the iD3Q12 MRT and D3Q13 MRT models under different relaxation time τ. The maximal pressure drop of the channel is increased and ${\rm{\text{δ}}}{x}={1}/{20}$ is fixed. The blank indicates that the computation is divergent.

2 -->

5.3.三维顶盖驱动的方腔流

三维方腔流包含旋涡运动, Ku等^[41]采用伪普方法计算的结果常被用作新数值方法模拟方腔流精度的评判标准, 因此我们也采用Ku的计算结果来检验iD3Q12 MRT模型的精度.
流动是在一个立方体盒子中进行的, 如图6所示. 流体是由最顶端的盖子以常速度$ U_0 = 1.0 $移动而被驱动. 流动满足三维不可压N-S方程. 雷诺数可由$ Re = U_0 L/\nu $计算, 这里$ L = 1.0 $是立方体盒子的长度, ν是剪切黏度. 在模拟中, 初始状态的速度和压力都设置为0. 在$ z = 0.5 $的截面上, 对称的边界条件设置为
图 6 三维顶盖驱动的方腔流示意图
Figure6. The schematic of three-dimensional lid-driven cavity flow

$ \frac{{\partial u_{x}}}{{\partial z}} = \frac{{\partial u_{y}}}{{\partial z}} = \frac{{\partial p}}{{\partial z}} = 0, \;{u_{z}} = 0. $

在$ y = 1 $的截面上$ x = 0 $和$ x = 1 $处$ {{u}} = 0 $. 在$ y = 1 $且$ x = 0 $和$ x = 1 $旁边的点处分别设置$ u_{x} = 0.3 $和$ u_{x} = 1.0 $. 初始条件和边界条件的设置与Ku等^[41]的相同. 此外, 为了满足低Ma假设, 固定$ c = 10 $, $ c \neq 1 $时MRT模型的计算见参考文献[31]. 速度域的收敛准则是:

$\begin{split} & {{\displaystyle \sum\nolimits_i \displaystyle\sum\nolimits_k {| {{{{u}}_{k}}({ x_i},t + 1000\text{δ} t) - {{{u}}_{k}}({ x_i},t)} |} } \over { \displaystyle\sum\nolimits_i \sum\nolimits_k {| {{{{u}}_{k}}({ x_i},t + 1000{\rm{\text{δ}}} t)} |} }} \\ & \leqslant 1.0\times10^{-14}, \end{split} $

对i和k的求和遍及所有的网格点和所有的速度方向. 首先采用iD3Q12 MRT模型模拟了Re = 100, 400和1000时的流动, 并将模拟结果与已有的Ku等^[41]的结果作对比, 如图7所示. 可以看出由iD3Q12 MRT模型计算的数值结果与Ku等^[41]的结果符合得很好, 因此我们提出的iD3Q12 MRT模型在模拟三维顶盖驱动的方腔流时是准确的.
图 7 不同的雷诺数下模拟方腔流, 在截面$z=0.5$处竖直和水平中心线的速度分布　(a) Re = 100; (b) $Re=400$; (c) $Re=1000$
Figure7. The velocity distribution in the vertical and horizontal center lines at section $z=0.5$ for cavity flows at different $Re$: (a) $Re=100$; (b) $Re=400$; (c) $Re=1000$.

不断地增大雷诺数时我们观察并记录了iD3Q12 MRT和He-Luo D3Q13 MRT模型在模拟方腔流时的敛散性, 结果显示在表6中. 从中可以看到在增大雷诺数时iD3Q12 MRT模型的稳定性比He-Luo D3Q13 MRT模型的稳定性稍弱一些.

Re	Model
	iD3Q12 MRT	He-Luo D3Q13 MRT
100	$\checkmark$	$\checkmark$
400	$\checkmark$	$\checkmark$
1000	$\checkmark$	$\checkmark$
1500	$\checkmark$	$\checkmark$
1600	$\checkmark$	$\checkmark$
1700	divergent	$\checkmark$
1800	divergent	divergent

表6不断增大雷诺数比较iD3Q12 MRT和He-Luo D3Q13 MRT模型在模拟方腔流时的稳定性. $\checkmark$代表收敛, 收敛准则是(39)式
Table6.Comparing the stability of iD3Q12 MRT and He-Luo D3Q13 MRT models for three-dimensional cavity flows when the Reynolds number is continuously increased. The tick represents convergence, the convergence criterion is formula (39).

提出了求解三维不可压缩流的12速多松弛格子Boltzmann模型, 即iD3Q12 MRT模型, 它可能是现有的能推导出不可压N-S方程的三维MRT模型中离散速度方向最少的一个模型, 因此原则上具有更高的计算效率. 构建iD3Q12 MRT模型的基本思想是构造一个12 × 12阶的正交变换矩阵, 并选取合适的矩平衡态使新模型能通过C-E展开恢复到不可压的N-S方程. 用iD3Q12 MRT模型分别模拟了稳态的泊肃叶流、非稳态的脉动流、顶盖驱动的方腔流, 验证了模型的准确性和稳定性, 并将其准确性和稳定性与He-Luo D3Q13 MRT模型作了对比.
对于泊肃叶流和脉动流, iD3Q12 MRT模型的数值解和解析解符合得很好, 这说明iD3Q12 MRT模型在模拟稳态流和非稳态流时都是准确的. 在不同的空间步长$ \Delta{x} $和松弛因子$ \lambda_{\nu} $下用iD3Q12 MRT模型和He-Luo D3Q13 MRT模型模拟泊肃叶流, 发现两个模型的全局相对误差$ {\rm GRE}_u $完全相同. 用两个模型模拟了脉动流并计算了不同空间步长$ \Delta{x} $和不同时刻下的全局相对误差, 发现iD3Q12 MRT模型的全局相对误差$ {\rm GRE}_u $比He-Luo D3Q13 MRT模型稍小一些. 全局相对误差随空间步长的变化表明两个模型在模拟稳态的泊肃叶流和非稳态的脉动流时都具有二阶精度. 还通过增大管道最大压降$ \Delta{p} $并调节空间步长$ \Delta{x} $或松弛时间τ的方式模拟了脉动流, 发现在大部分的参数下iD3Q12 MRT模型计算的流场的全局相对误差$ {\rm GRE}_u $比He-Luo D3Q13 MRT模型计算的流场的全局相对误差要小一些, 但随着最大压降的增大iD3Q12 MRT模型比D3Q13 MRT模型先发散. 这说明在模拟非稳态流时iD3Q12 MRT模型比D3Q13 MRT模型的准确性稍高但在稳定性方面弱一些. 此外, 用iD3Q12 MRT模型模拟了包含旋涡运动的三维顶盖驱动方腔流, 发现iD3Q12 MRT模型的模拟结果与已有的Ku等的结果符合得很好, 并且iD3Q12 MRT模型能模拟的最大雷诺数比He-Luo D3Q13 MRT模型的要稍小.

通过C-E展开将iD3Q12 MRT模型恢复到三维不可压的N-S方程. 对于不可压缩流有

$\small O({\rm{\text{δ}}} p) = O({\rm{\text{δ}}} \rho ) = O({Ma^2}),\tag{A1a}$

$\small O( u) = O(Ma),\tag{A1b}$

这里Ma代表马赫数, $ {\rm{\text{δ}}} p $和$ {\rm{\text{δ}}} \rho $分别代表压力和密度的脉动量.
首先引入展开

$ {f_i }( x + {{ c}_i }{\text{δ}t},t + {\text{δ}t}) = \sum\limits_{n = 0}^\infty {\frac{{{\varepsilon ^n}}}{{n!}}} {D_{t}^n}{f_i }( x,t),\tag{A2a}$

$ {f_i } = \sum\limits_{n = 0}^\infty {{\varepsilon ^n}f_i^{(n)}},\tag{A2b}$

$ {\partial _t} = \sum\limits_{n = 0}^\infty {{\varepsilon ^n}{\partial _{{t_n}}}},\tag{A2c}$

这里$ \varepsilon = {\rm{\text{δ}}} t $, $ {D_i} \equiv {\partial _t} \!+\! {{ c}_i } \cdot \nabla = {\partial _t} \!+\! {{ c}_{{i}{\alpha}} } \cdot {\partial_{\alpha}} = {\partial _t}\!+\! c_{i x}\partial_x+$$c_{i y}\partial_y+c_{i z}\partial_z $, 运用上述(A2)式展开, 并将演化方程

$\begin{split} &{f_i }({{x}} + {{{c}}_i }{\text{δ}t},t + {{\text{δ}} t}) - {f_i }({{x}},t)\\ =\, & - \sum\limits_{j = 1}^{12} {{\varLambda _{i j}}({f_j}({{x}},t) - f_j^{({\rm{eq}})}({{x}},t))} , \;\;i =1\text{---}12, \end{split}\tag{A3}$

做Taylor展开, 按照ε的不同阶可得:

$ O({\varepsilon ^0}): f_i ^{(0)} = f_i ^{({\rm{eq}})},\tag{A4a}$

$ O({\varepsilon ^1}): {D_{t_0}}f_i ^{(0)} = - \sum\limits_{j = 1}^{12} {{{\varLambda}_{i j}}f_j^{(1)}},\tag{A4b}$

$ O({\varepsilon ^2}): {\partial _{{t_1}}}f_i ^{(0)} + {D_{t_0 }}\left( {f_i ^{(1)}\! -\! \frac{1}{2}\sum\limits_{j = 1}^{12} {{\varLambda _{i j}}f_j^{(1)}} } \right) \!=\! - \sum\limits_{j = 1}^{12} {{{\varLambda}_{i j}}f_j^{(2)}}, \tag{A4c}$

这里$ {D_{i_0}} \equiv {\partial _{{t_0}}} + {{ c}_i } \cdot {\nabla} = {\partial _t} + {{ c}_{{i}{\alpha}} } \cdot {\partial_{\alpha}} = {\partial _{t_0}}+c_{i x}\partial_{x}+ $$c_{i y}\partial_{y}+c_{i z}\partial_{z} $, 并且将方程(A4b)代入到方程(A4c)中, 方程(A4)可以转化到矩空间中:

$ {{ m}^{(0)}} = {{ m}^{({\rm{eq}})}}, \tag{A5a}$

$ ({\partial _{{t_0}}} { I} + {\hat {{ C}_k}}{\partial _k}){{ m}^{(0)}} = - \hat {\boldsymbol {\varLambda} }{{ m}^{(1)}},\tag{A5b}$

$ {\partial _{{t_1}}}{{ m}^{(0)}} + ({\partial _{{t_0}}} { I} + {\hat {{ C}_k}}{\partial _k})\left( { I - \frac{1}{2} \hat{ \boldsymbol {\varLambda} }} \right){{ m}^{(1)}} = - \hat{\boldsymbol {\varLambda}} {{ m}^{(2)}},\tag{A5c}$

$ {\hat {{ C}_k}} = { T}{{ C}_k}{{ T}^{ - 1}} $, $ { C}_k $是以离散速度$ { c}_i $的第k个分量作为对角线元素的对角矩阵,

$ \begin{split}{{ m}^{(n)}} =\, & \left(0,0,0,0,3S_{xx}^{(n)},S_{{\omega}{\omega}}^{(n)}, S_{xy}^{(n)}, S_{yz}^{(n)},\right.\\ & \left.S_{zx}^{(n)},h_x^{(n)},h_y^{(n)}, h_z^{(n)}\right)^{\prime}, \; \; n = 1,2. \end{split}\tag{A6}$

值得注意的是, $ P_1 $和$ {u_{x}}, {u_{y}}, {u_{z}} $在低Ma下是守恒量, 这样${{P}_{1}}^{(n)} $和$ u_{x}^{(n)}, u_{y}^{(n)}, u_{z}^{(n)} $在表达式(A6)中都为0. 展开方程(A5b)有

$\begin{split}&{\partial _{{t_0}}}\left[\!\! {\begin{array}{*{20}{c}} {u^2/2 + 3p/2}\\ {{u_x}}\\ {{u_y}}\\ {{u_z}}\\ {3u_x^2-u^2}\\ {u_y^2-u_z^2}\\ {{u_x}{u_y}}\\ {{u_y}{u_z}}\\ {{u_x}{u_z}}\\ 0\\ 0\\ 0 \end{array}} \!\!\right] + {\partial _{x}}\left[\!\! {\begin{array}{*{20}{c}} {{u_x}}\\ {u_x^2+p}\\ {{u_x}{u_y}}\\ {{u_x}{u_z}}\\ {{u_x}}\\ 0\\ {{u_y}/2}\\ 0\\ {{u_z}/2}\\ {u_y^2-u_z^2}\\ {-{u_x}{u_y}}\\ {{u_x}{u_z}} \end{array}} \!\!\right] + {\partial _{y}}\left[\!\! {\begin{array}{*{20}{c}} {{u_y}}\\ {{u_x}{u_y}}\\ {u_y^2+p}\\ {{u_y}{u_z}}\\ {-{u_y}/2}\\ {{u_y}/2}\\ {{u_x}/2}\\ {{u_z}/2}\\ 0\\ {{u_x}{u_y}}\\ {u_z^2-u_x^2}\\ {-{u_y}{u_z}} \end{array}} \!\!\right]\\ +\,& {\partial _{z}}\left[\!\! {\begin{array}{*{20}{c}} {{u_z}}\\ {{u_x}{u_z}}\\ {{u_y}{u_z}}\\ {u_z^2+p}\\ {-{u_z}/2}\\ {-{u_z}/2}\\ 0\\ { {u_y}/2}\\ { {u_x}/2}\\ {-{u_x}{u_z}}\\ {{u_y}{u_z}}\\ {u_x^2-u_y^2} \end{array}} \!\!\right] = - \left[\!\! {\begin{array}{*{20}{c}} 0\\ 0\\ 0\\ 0\\ {-{\lambda_{\nu} }3S_{xx}^{(1)}}\\ {-{\lambda_{\nu} }S_{\omega \omega }^{(1)}}\\ {-{\lambda_{\nu}^{\prime} }S_{xy}^{(1)}}\\ {-{\lambda_{\nu}^{\prime} }S_{yz}^{(1)}}\\ {-{\lambda_{\nu}' }S_{zx}^{(1)}}\\ {-{\lambda_t}h_{x}^{(1)}}\\ {-{\lambda_t}h_{y}^{(1)}}\\ {-{\lambda_t}h_{z}^{(1)}} \end{array}} \!\!\right]\end{split}\tag{A7}.$

(A7)式的前4个方程是

$ {\partial _{{t_0}}}(3p/2 + { u}^2/2)+{\partial _{x}}{u_x} + {\partial _{y}}{u_y} + {\partial _{z}}{u_z} = 0,\tag{A8}$

$ {\partial _{{t_0}}}{u_x} + {\partial _{x}}(p + u_x^2) + {\partial _{y}}({u_x}{u_y}) + {\partial _{z}}({u_x}{u_z}) = 0,\tag{A9a}$

$ {\partial _{{t_0}}}{u_y} + {\partial _{x}}({u_x}{u_y}) + {\partial _{y}}(p + u_y^2) + {\partial _{z}}({u_y}{u_z}) = 0,\tag{A9b}$

$ {\partial _{{t_0}}}{u_z} + {\partial _{x}}({u_x}{u_z}) + {\partial _{y}}({u_y}{u_z}) + {\partial _{z}}(p + u_z^2) = 0.\tag{A9c}$

从(A1a)与(A1b) 可以得到$ O({\rm{\text{δ}}} p) = O(M^2) $和$ O( u) = $ O(M), 有$ {\partial _{{t_0}}}(3 p/2 + u^2/2) = O(M^2) $. 忽略$ O(M^2) $项, 方程(A8)变成

$ {\partial _{x}}{u_x} + {\partial _{y}}{u_y} + {\partial _{z}}{u_z} = 0,\tag{A10}$

它是不可压N-S方程的连续方程. 从方程(A5c)中可得

$ \begin{split} {\partial _{{t_1}}}{u_x} \, &+ {\partial _{x}}[(1 - {\lambda_{\nu}}/2)S_{xx}^{(1)}] + {\partial _{y}}[(1 - {\lambda_{\nu}^{\prime} }/2)S_{xy}^{(1)}] \\ &+ {\partial _{z}}[(1 - {\lambda_{\nu}^{\prime} }/2)S_{zx}^{(1)}] = 0, \end{split}\tag{A11a}$

$ \begin{split} {\partial _{{t_1}}}{u_y}\, & + {\partial _{x}}[(1 - {\lambda_{\nu}^{\prime} }/2)S_{xy}^{(1)}] + {\partial _{y}}[(1 - {\lambda_{\nu} }/2)(S_{\omega \omega }^{(1)} \\&- S_{xx}^{(1)})/2] + {\partial _{z}}[(1 - {\lambda_{\nu}^{\prime} }/2)S_{yz}^{(1)}] = 0, \end{split}\tag{A11b} $

$ \begin{split} {\partial _{{t_1}}}{u_z}\, & +{\partial_{x}}[1-{\lambda_{\nu}^{\prime}}/2)S_{zx}^{(1)}]-{\partial_{z}}[(1-{\lambda_{\nu}}/2)(S_{\omega\omega}^{(1)}\\ &+S_{xx}^{(1)})/2] + {\partial _{y}}[(1 - {\lambda_{\nu}^{\prime} }/2)S_{yz}^{(1)}] = 0, \end{split}\tag{A11c}$

由(A9)式 + $ \varepsilon\times $(A11)式可得

$ \begin{split} &{\partial _t}{u_x} + {\partial _x}(p + u_x^2) + {\partial _y}({u_x}{u_y}) + {\partial _z}({u_x}{u_z}) \\=\,& - \varepsilon \{{\partial _x}[(1 - {\lambda_{\nu}}/2)S_{xx}^{(1)}]\\ &+{\partial _y}[(1 - {\lambda_{\nu}^{\prime} }/2)S_{xy}^{(1)}]+{\partial_z}[(1-{\lambda_{\nu}^{\prime}}/2)S_{zx}^{(1)}]\}, \end{split} \tag{A12a}$

$ \begin{split} & {\partial _t}{u_y} + {\partial _x}({u_x}{u_y}) + {\partial _y}(p + u_y^2) + {\partial _z}({u_y}{u_z}) \\=\, & - \varepsilon\{ {\partial_x}[(1 - {\lambda_{\nu}^{\prime}}/2)S_{xy}^{(1)}]\\&+ {\partial_y}[(1-{\lambda_{\nu}}/2)(S_{\omega \omega}^{(1)}-S_{xx}^{(1)})/2]+ {\partial_z}{[(1 - {\lambda_{\nu}^{\prime} }/2)S_{yz}^{1}]}\}, \end{split}\tag{A12b}$

$ \begin{split} &{\partial _t}{u_z} + {\partial _x}({u_x}{u_z}) + {\partial _y}({u_y}{u_z}) + {\partial _z}(p + u_z^2)\\ =\, & - \varepsilon\{ {\partial _x}[(1 - {\lambda_{\nu}^{\prime} }/2)S_{zx}^{(1)}]\\ & +{\partial _y}[(1 - {\lambda_{\nu}^{\prime} }/2)S_{yz}^{(1)}] -{\partial _z}[ (1 - {\lambda_{\nu}}/2)(S_{\omega \omega}^{(1)}+S_{xx}^{(1)})/2]\}, \end{split}\tag{A12c}$

根据方程(A7), 算出

$ S_{xx}^{(1)} = - \frac{1}{3{\lambda}_{\nu}} \left[ {{\partial_{{t_0}}}(3u_x^2-u^2)+{\partial_{x}}{u_x}- \frac{1}{2}{\partial_{y}}{u_y}-\frac{1}{2}{\partial_{z}{u_z}}} \right],\tag{A13a}$

$ S_{\omega \omega }^{(1)} = - \frac{1}{{\lambda}_{\nu}}\left[ {{\partial_{{t_0}}}(u_y^2-u_z^2) +\frac{1}{2}{\partial_{y}}{u_y}-\frac{1}{2}{\partial_{z}}{u_z}} \right], \tag{A13b} $

$ S_{xy}^{(1)} = - \frac{1}{{\lambda}_{\nu}^{\prime} } \left[ {{\partial _{{t_0}}}({u_x}{u_y}) + \frac{1}{2}{\partial _{x}}{u_y} + \frac{1}{2}{\partial _{y}}{u_x}} \right], \tag{A13c}$

$ S_{yz}^{(1)} = - \frac{1}{{\lambda}_{\nu}^{\prime}}\left[ {{\partial _{{t_0}}}({u_y}{u_z}) + \frac{1}{2}{\partial _{y}}{u_z} + \frac{1}{2}{\partial_{z}{u_y}}} \right], \tag{A13d}$

$ S_{zx}^{(1)} = - \frac{1}{{\lambda}_{\nu}^{\prime} }\left[ {{\partial _{{t_0}}}({u_x}{u_z}) + \frac{1}{2}{\partial _{x}}{u_z} + \frac{1}{2}{\partial _{z}}{u_x}} \right]. \tag{A13e}$

使用方程(A9), 得到

$ {\partial _{{t_0}}}{u_x^2} = -2u_x \big[{\partial _{x}}(p + u_x^2) + {\partial _{y}}({u_x}{u_y}) + {\partial _{z}}({u_x}{u_z}) \big], \tag{A14a} $

$ {\partial _{{t_0}}}{u_y^2} = -2u_y \big[{\partial _{x}}({u_x}{u_y}) + {\partial _{y}}(p + u_y^2) + {\partial _{z}}({u_y}{u_z})\big],\tag{A14b}$

$ {\partial _{{t_0}}}{u_z^2} = -2u_z \big[{\partial _{x}}({u_x}{u_z}) + {\partial _{y}}({u_y}{u_z}) + {\partial _{z}}(p + u_z^2)\big].\tag{A14c}$

从方程(A1)可知$ \partial_{t_0}u_x^2 $, $ \partial_{t_0}u_y^2 $和$ \partial_{t_0}u_z^2 $都是$ O(M^3) $项. 忽略$ O(M^3) $项和$ \partial_{t_0}u^2 $项, 方程(A13)可以写作

$ S_{xx}^{(1)} = - \frac{1}{3{\lambda}_{\nu}}\left( {{\partial_{x}}{u_x} -\frac{1}{2}{\partial_{y}}{u_y}-\frac{1}{2}{\partial_{z}}{u_z}} \right),\tag{A15a}$

$ S_{\omega \omega }^{(1)} = - \frac{1}{2{\lambda}_{\nu} }({\partial _{y}}{u_y} - {\partial _{z}}{u_z}), \tag{A15b}$

$ S_{xy}^{(1)} = - \frac{1}{2{\lambda}_{\nu}^{\prime} }({\partial _{x}}{u_y} + {\partial _{y}}{u_x}),\tag{A15c}$

$ S_{yz}^{(1)} = - \frac{1}{2{\lambda}_{\nu}^{\prime}}({\partial _{y}}{u_z} + {\partial_{z}}{u_y}), \tag{A15d}$

$ S_{zx}^{(1)} = - \frac{1}{2{\lambda}_{\nu}^{\prime} }({\partial _{x}}{u_z} + {\partial _{z}}{u_x}),\tag{A15e}$

将方程(A15)代入(A12)式, 可得

$ \begin{split} & {\partial _t}{u_x} + {\partial _x}(u_x^2) + {\partial _y}({u_x}{u_y}) + {\partial _z}({u_x}{u_z})\\ =\, & -{\partial _x}p+\nu (\partial _x^2{u_x} + \partial _y^2{u_x} + \partial _z^2{u_x})\\ &+\xi ({\partial_x}({\partial_x}(u_x)+{\partial_y}(u_y)+{\partial_z}(u_z))), \end{split}\tag{A16a} $

$ \begin{split}& {\partial _t}{u_y} + {\partial _x}({u_x}{u_y}) + {\partial _y}(u_y^2) + {\partial _z}({u_y}{u_z})\\ =\, & -{\partial _y}p+\nu (\partial _x^2{u_y} + \partial _y^2{u_y} + \partial _z^2{u_y})\\ & +\xi({\partial_y}({\partial_x}(u_x)+{\partial_y}(u_y)+{\partial_z}(u_z))), \end{split}\tag{A16b}$

$ \begin{split} & {\partial _t}{u_z} + {\partial _x}({u_x}{u_z}) + {\partial _y}({u_y}{u_z}) + {\partial _z}(u_z^2) \\ =\, & -{\partial _z}p+\nu (\partial _x^2{u_z} + \partial _y^2{u_z}+ \partial _z^2{u_z})\\ & +\xi({\partial_z}({\partial_x}(u_x)+{\partial_y}(u_y)+{\partial_z}(u_z))), \end{split}\tag{A16c}$

其中,

$\begin{split} & \nu = \dfrac{1}{4} \left( {\tau- \dfrac{1}{2}} \right){{\rm{\text{δ}}} t} = \dfrac{1}{2} \left( {\dfrac{1}{{\lambda}_{\nu}^{\prime}}-\dfrac{1}{2}} \right){{\rm{\text{δ}}} t}, \\ & \xi = \dfrac{1}{12}\left( {{\tau-\dfrac{1}{2}}} \right){\text{δ} t}, ~\tau = \dfrac{1}{\lambda_{\nu} } \end{split} $

由(A10)式, 省略方程(A16)中$ {\xi}(\nabla({\nabla}\cdot{{{u}}})) $和$ { u}{\nabla}\cdot{ u} $两项, 得

$ \begin{split} & {\partial _t}{u_x} + {u_x}{\partial _x}{u_x} + {u_y}{\partial _y}{u_x} + {u_z}{\partial _z}{u_x} \\=\, & -{\partial _x}p+\nu (\partial _x^2{u_x} + \partial _y^2{u_x} + \partial _z^2{u_x}), \end{split}\tag{A17a}$

$ \begin{split} & {\partial _t}{u_y} + {u_x}{\partial _x}{u_y} +{u_y}{\partial _y}{u_y} + {u_z}{\partial _z}{u_y} \\ = \, & -{\partial _y}p+\nu (\partial _x^2{u_y} + \partial _y^2{u_y} + \partial _z^2{u_y}), \end{split}\tag{A17b} $

$ \begin{split} &{\partial _t}{u_z} + {u_x}{\partial _x}{u_z} + {u_y}{\partial _y}{u_z} + {u_z}{\partial _z}{u_z} \\=\, & -{\partial _z}p+\nu (\partial _x^2{u_z} + \partial _y^2{u_z} + \partial _z^2{u_z}), \end{split}\tag{A17c}$

这就是不可压N-S方程动量方程的分量形式. 至此, 在低马赫数假设下, 已经通过C-E展开将iD3Q12 MRT模型恢复到不可压的N-S方程, 它可以写作向量的形式, 见方程(25).