LOCAL GRID REFINEMENT APPROACH FOR LATTICE BOLTZMANN METHOD: DISTRIBUTION FUNCTION CONVERSION BETWEEN COARSE AND FINE GRIDS
-
摘要: 格子Boltzmann方法作为一种高效的介观计算流体力学方法在过去20多年里得到快速发展, 其相对较高的计算效率和灵活性使其可以适用于各种复杂流动的模拟. 然而标准的格子Boltzmann方法只能使用均匀的直角网格, 这种网格排布方式并不利于复杂流动的计算. 为此, 基于格子Boltzmann方法的局部网格加密算法在文献中被提出. 该算法需要在局部加密的界面处将粗细网格间的分布函数转换后交换. 目前分布函数的转换方式大多是在没有源项的情况下推导的, 而且现存考虑源项时转换公式的推导也都是基于Chapman-Enskog展开; 其推导过程相对复杂, 且需要对分布函数的非平衡态部分做一阶Chapman-Enskog近似, 这有可能会限制局部网格加密算法在高阶格子Boltzmann方法中的应用. 文章在忽略时空离散误差的前提下, 以保证连续分布函数变量以及物理松弛系数一致为基础, 构建了一套规范且简洁的粗细网格间在考虑任意源项时, 分布函数转换关系的推导过程, 该方法不依赖于Chapman-Enskog展开以及Chapman-Enskog近似, 且该方法既可以适用于单松弛碰撞模型也可以适用于多松弛碰撞模型. 此外, 还从理论上证明了, 保证粗细网格间非平衡态部分的一阶 Chapman-Enskog 近似一致, 便可以保证整个非平衡态部分的一致, 这将有助于扩展局部网格加密算法中转换关系的应用范围. 最后, 通过对强迫泰勒−格林涡流动、平板泊肃叶流中对流−扩散问题和顶盖驱动方腔流动进行数值模拟, 良好的数值结果证实了转换关系对复杂源项的适应性以及局部网格加密技术在处理复杂流动问题方面的优势. 同时, 通过对一维剪切波问题的模拟, 发现由局部网格加密引起的数值黏性与加密区域的选取有很大的关系.
-
关键词:
- 格子Boltzmann方法 /
- 局部网格加密 /
- 数据转换
Abstract: Lattice Boltzmann method, as an efficient mesoscopic computational fluid dynamics method, has been developed rapidly in the past two decades. Its relatively high computational efficiency and flexibility make it suitable for the simulation of various complex flows. However, due to its own limitations, the standard lattice Boltzmann method typically utilize uniform rectangular grid, which is not suitable for the simulation of complex flows. Therefore, local grid refinement based on the lattice Boltzmann method has been considered by many researchers. For this purpose, the distribution functions between coarse and fine grids need to be converted at the interface of coarse and fine grids. At present, most conversion methods of distribution function are derived without the presence of the source term, and the previously limited derivation of conversion formulas considering the source term was based on the first-order Chapman-Enskog expansion, which is relatively complex and may limit the application of local grid refinement algorithm in higher-order lattice Boltzmann methods. In this paper, we provide a concise derivation to relate the distribution functions between coarse and fine grids considering an arbitrary source term, based on consistency requirements of the distribution function of the continuous Boltzmann equation between coarse and fine grids. The proposed method is independent of the Chapman-Enskog expansion and Chapman-Enskog approximation, and can be applied to both single relaxation time and multiple relaxation time collision models. In addition, this paper also proves theoretically that the consistency of the first-order Chapman-Enskog approximation of the non-equilibrium distribution between the coarse and fine girds can ensure the consistency of the entire non-equilibrium distribution, which expands the applicability of the previous conversion relationship. Finally, these theoretical results are validated by numerical simulations of a forced Taylor-Green vortex flow, convection-diffusion in a planar Poiseuille flow and lid-driven cavity flow. The good numerical results confirm the adaptability of the conversion relation in the presence of complex source terms and the advantages of local grid refinement technology in dealing with complex flow problems. At the same time, through the simulation of one-dimensional shear wave problem, it is found that the numerical viscosity caused by local grid refinement has a great relationship with the selection of refinement region.-
Keywords:
- lattice Boltzmann method /
- grid refinement /
- data conversion
-
引言
在实际的工程应用中, 各种各样的流体力学问题随处可见, 而如何准确且高效的描述其对应的流动特征一直是工业界和科学界重点关注的问题. 近半个世纪以来, 计算机技术的迅猛发展使得大规模数值求解复杂问题成为可能, 而数值模拟作为一种有效的研究手段也得到了更加广泛的关注和应用. 不同于直接数值求解纳维−斯托克斯(Navier-Stokes, N-S)方程的传统计算流体力学方法, 介观计算流体动力学方法是直接数值求解基于分子动力学理论的Boltzmann方程, 而格子玻尔兹曼方法(lattice Boltzmann method, LBM)则是近30年来出现的一种较为高效的介观计算流体动力学方法. 它较高的计算效率和灵活性使其可以适用于各种高度复杂的物理现象, 如湍流[1]、燃烧[2]、多相流[3]、磁流体[4]和粒子悬浮[5-6]等问题. 由于其数值耗散较低, 故LBM也常被应用于气动声学模拟中[7-8]. 同时, LBM也凭借演化过程中“碰撞−迁移”的局部性, 使得其非常适合于大规模的并行计算[9-10]. 除此之外, 许多学者为了充分利用LBM的优点, 还通过反设计的思想在原始的Boltzmann方程中添加一些特定的源项, 以使得修改后的LBM可以满足特定需求[11-12], 而这也使得LBM方法的适用范围得到了极大的扩展.
流体的流动特点一般是随着Reynolds数的增加, 会产生越来越精细的动力学结构. 而工业领域里涉及到的流动问题大多是高Reynolds数的湍流, 通过直接数值计算来解析所有的湍流结构往往需要较高的计算成本. 因此, 为了在有限的计算资源下获得较为准确的结果, 往往需要提高局部分辨率来解析这些局部流域内包含的精细流动结构. 但由于LBM算法在时空离散时是沿着特征线进行的离散, 故在给定离散速度集后LBM中的时间步长和空间步长便通过离散速度相互耦合, 这一特点导致标准的LBM只能使用均匀的直角网格.
为了克服标准LBM只能使用均匀网格这一缺陷, Filippova等[13]提出了基于LBM的局部网格加密(FH)方法. 该方法首先需要使用均匀的直角粗网格(背景网格)覆盖整个流场, 然后在需要进行加密的区域上对网格进行加密, 不同网格分辨率区域下仍使用标准的LBM演化关系且使用同一离散速度集. 而粗细网格之间的耦合则通过流场物性参数以及分布函数间的转换实现, 如可以通过保证粗细网格间运动黏度系数一致来得到粗细网格间格子松弛系数的关系, 然后通过保证粗细网格间宏观量一致得到分布函数的转换关系. 文中以二维圆柱扰流为测试算例, 在使用局部网格加密技术后得到了较好的数值结果, 这也验证了网格局部加密的有效性. 局部网格加密算法的提出对于标准LBM的意义是重大的, 这意味着其可以通过局部加密来节省大量的计算成本, 也使得LBM被大规模应用于工业界成为了可能.
与FH方法不同, Lin等[14]建议不对粗细网格间分布函数的非平衡部分做任何缩放, 而是直接交换粗细网格的分布函数. Dupuis等[15]比较了FH方法和Lin的方法后发现, 当考虑非平衡部分的转换时, 可以得到更精确的数值结果. 而Touil等[16]更是直接指出, Lin的方法实际上是混淆了连续分布函数和离散分布函数之间的区别, 其直接在粗细网格间交换分布函数会导致该方法的精度下降. 此外, Rohde等[17]提出了一种基于有限体积思想且可以适用于任何碰撞模型的局部网格加密算法, 但其也犯了与Lin相同的错误, 即不对分布函数做任何转换, 这种处理方式在模拟高Reynolds数下的流动时, 很容易在网格加密界面处出现间断和非物理震荡. 而许多学者[18-21]的相关研究也表明, 在局部网格加密算法中, 粗细网格间非平衡态分布函数的转换是必要且不可或缺的.
众所周知, LBM的灵活性很大程度上来自于可以通过Chapman-Enskog (CE)展开[22]灵活设计源项. 因此, 在考虑源项的前提下, 建立粗细网格间分布函数的转换关系, 对于扩展局部网格加密算法在LBM中的应用范围至关重要. Huang等[23]基于CE展开, 提出了一种在多松弛碰撞模型(multiple relaxation time, MRT)[24]下利用多层网格技术计算被动标量的模型. 该方法在考虑外力项存在的情况下, 在矩空间下执行粗细网格间的数据转换. Liou等[25]同样也基于CE展开, 提出了一种在BGK (Bhatnagar-Gross-Krook)模型[26]下, 同时考虑外力和标量源项的加密方法. 值得注意的是, 这些方法都是在 CE 展开的基础上导出的, 其推导过程较为复杂, 且在推导中对分布函数的非平衡态使用了一阶CE近似, 而这有可能限制局部网格加密算法在高阶LBM [27]中的应用. 本文的主要目标是, 在考虑任意源项存在的前提下, 提出一套规范且简洁的粗细网格在重合点处的数据转换关系推导方法.
本文的其余部分按如下顺序展开. 在第 1 节, 梳理BGK模型和MRT模型下LBM的时空离散过程, 尤其是指出在时空离散过程中离散分布函数变量的引入, 从而为后文局部网格加密算法中的数据转换进行必要的铺垫; 同时, 梳理BGK模型和MRT模型下非平衡态分布函数的初始化方法, 为后文数值验证中分布函数的初始化提供依据. 第 2 节则是凭借第 1 节中的理论基础, 分别详细给出BGK模型和MRT模型下粗细网格间数据转换关系的推导过程. 在第 3 节中, 通过对强迫泰勒−格林涡流动和平板泊肃叶流中对流−扩散问题进行数值模拟, 以验证第 2 节中转换关系对复杂源项的适应性; 通过对顶盖驱动方腔流动进行数值模拟来展示局部网格加密技术在处理复杂流动问题方面的优势; 通过对一维剪切波问题的模拟, 来观察局部网格加密技术对LBM在数值耗散方面的影响. 最后, 在第 4 节给出本文结论.
1. 理论基础
速度离散后的Boltzmann 方程(discrete velocity Boltzmann equation, DVBE)可以表示如下[28]
$$ \partial_t f_\alpha({\boldsymbol{x}},t)+{\boldsymbol{\xi}}_\alpha \cdot {\boldsymbol{\nabla}} f_\alpha({\boldsymbol{x}},t)=\varOmega_\alpha({\boldsymbol{x}},t)+S_\alpha({\boldsymbol{x}},t) $$ (1) 这里, $ {\boldsymbol{\xi}}_\alpha $ 表示第 $ \alpha $ 个离散速度点对应的离散速度; $ f_\alpha({\boldsymbol{x}},t) $代表$ t $时刻在位置$ {\boldsymbol{x}} $处以速度$ {\boldsymbol{\xi}}_\alpha $移动的粒子对应的分布函数; $\varOmega_{\alpha}$代表碰撞项; $ S_{\alpha} $为可以代表包括体力项在内的任意源项. 以二维问题中被广泛使用的D2Q9速度集为例, 如图1 所示, 其离散速度点可以表示如下 [28]
$$ \begin{split} &\{\xi_\alpha\}_{\alpha=1}^{9}=\\ &\quad \sqrt{3}\left\{\begin{array}{ccccccccc} 0,&1,&0,&-1,&0,&1,&-1,&-1,&1 \\ 0,&0,&1,&0,&-1,&1,&1,&-1,&-1 \end{array}\right\}\end{split} $$ (2) 不同离散速度点对应的权重为
$$ w_\alpha=\left\{\begin{split} &4/9,\quad \alpha=1\\ &1/9, \quad \alpha=2,3,4,5\\ &1/36, \quad\alpha=6,7,8,9 \end{split}\right.$$ (3) 需要注意的是, 本文使用的D2Q9离散速度集长度未经缩放, 即格子声速: $ C_s=1. $
在标准的LBM中, 时空离散是沿着特征线进行的离散, 这种离散方法导致了时间步长和空间步长通过离散的速度集进行了耦合. 这就使得整个计算域都必须使用均匀的直角网格, 以保证流体粒子在一个时间步长$ \Delta t $后可以精确地移动到邻近的网格节点上, 这也是LBM需要使用 on-lattice 速度集的原因. 而对于其他的介观方法, 如: 基于有限体积法对DVBE经行离散的DUGKS [29]等方法, 这些方法中时间步长和空间步长是解耦的, 故其可以使用非均匀网格, 且既可以使用on-lattice的速度集也可以使用off-lattice的速度集. 接下来, 我们将给出标准LBM中对DVBE的时空离散过程
直接对DVBE左右两侧沿特征线积分, 可得如下结果
$$ \begin{split} &f_{\alpha}\left({\boldsymbol{x}}+{\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t\right)-f_{\alpha}({\boldsymbol{x}}, t)=\\ &\qquad \frac{\Delta t}{2}\left[\left(\varOmega_{\alpha}^{\prime} + S_{\alpha}^{\prime}\right)+\left(\varOmega_{\alpha}({\boldsymbol{x}},t) + S_{\alpha}({\boldsymbol{x}},t)\right)\right] + O(\Delta t^3)\end{split} $$ (4) 从上式不难看出, 分布函数$ f_\alpha $从$ t $时刻到$ t+\Delta t $时刻, 时间上推进了$ \Delta t $而空间上推进的长度为$ \Delta {\boldsymbol{x}}={\boldsymbol{\xi}}_\alpha\Delta t $, 这也正是LBM中时空耦合的原因. 同时在推导上式的过程中, 我们对DVBE等号右边项使用了梯形积分公式来近似, 其中$ \varOmega_{\alpha}^{\prime},S_{\alpha}^{\prime} $分别代表$ \varOmega_{\alpha}\left({\boldsymbol{x}}+{\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t\right),S_\alpha\left({\boldsymbol{x}}+{\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t\right) $. 不难发现式 (4) 右边项中包含了$ t+\Delta t $时刻的未知量, 而引入新的变量$ g_\alpha({\boldsymbol{x}},t) $可以消除这种隐式 [30], 该离散分布函数变量的表达式如下
$$ g_\alpha({\boldsymbol{x}},t)=f_\alpha({\boldsymbol{x}},t) - \frac{\Delta t}{2}(\varOmega_{\alpha} + S_{\alpha}) $$ (5) 现在可以使用离散分布函数变量$ g_\alpha $对略去高阶小量的式 (4) 重新表示
$$ g_{\alpha}\left({\boldsymbol{x}}+{\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t\right)=g_{\alpha}({\boldsymbol{x}}, t)+\Delta t(\varOmega_{\alpha} + S_{\alpha}) $$ (6) 在给定碰撞项$ \varOmega_\alpha $以及源项$ S_\alpha $的具体形式后, 结合式 (5) 对式 (6) 简单整理后便可以得到对应碰撞模型下的LBM演化式. 下面以BGK碰撞模型和MRT碰撞模型为例, 给出两种模型下LBM演化式的详细推导过程.
1.1 BGK碰撞模型
速度离散后的BGK碰撞模型可以直接表示为
$$ \varOmega_{{\rm{BGK}},\alpha}({\boldsymbol{x}},t)=-\frac{1}{\tau}\left(f_\alpha({\boldsymbol{x}},t)-f^{eq}_\alpha({\boldsymbol{x}},t)\right) $$ (7) 式中, $ \tau $为物理松弛系数, 其与网格分辨率无关; 而$ f_\alpha^{eq} $为速度离散后的Maxwell 平衡态, 其对应的Hermite展开到前二阶的表达式(等温)如下[28]
$$ \begin{split} &f_\alpha^{eq}=w_\alpha \sum\limits_{n=0}^{2} \frac{1}{n!}{\boldsymbol{a}}^{eq,(n)}({\boldsymbol{x}},t){\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha)=\\ &\qquad w_\alpha\rho\left\{1+\underbrace{{\boldsymbol{\xi}}_\alpha\cdot{\boldsymbol{u}}}_{1{\rm{st}}}+\underbrace{\frac{1}{2}\left[({\boldsymbol{\xi}}_\alpha\cdot{\boldsymbol{u}})^2-u^2\right]}_{2{\rm{nd}} \ {\rm{order}}}\right\}\end{split} $$ (8) 式中, ${\boldsymbol{a}}^{eq,(n)}=\displaystyle\sum _{\alpha=1}^{q}f_\alpha^{eq} {\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha)$ 代表平衡态分布函数$ f^{eq}_\alpha $对应的第$ n $阶Hermite矩, 前两阶Hermite矩可以表示为
$$ a^{eq,(0)}=\rho,\qquad a^{eq,(1)}_{i}=\rho u_i,\qquad a^{eq,(2)}_{ij}=\rho u_i u_j $$ (9) 式中, $ q $ 为离散速度集中离散速度点的个数. 而$ {\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha) $代表以$ {\boldsymbol{\xi}}_\alpha $为自变量的第$ n $阶Hermite多项式[31], 前两阶的Hermite多项式可以表示如下
$$ \left.\begin{split} &H^{(0)}({\boldsymbol{\xi}}_\alpha)=1\\ &H^{(1)}_i({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,i}\\ &H^{(2)}_{ij}({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,i}\xi_{\alpha,j}-\delta_{ij} \end{split}\right\} $$ (10) 其中, 式 (8) 中的密度$ \rho $和速度$ {\boldsymbol{u}} $可以由$ f_\alpha $以及离散速度$ {\boldsymbol{\xi}}_\alpha $表示如下
$$ \rho=\sum\limits_{\alpha=1}^{q}f_\alpha,\quad \rho u_i=\sum\limits_{\alpha=1}^{q}f_\alpha \xi_{\alpha,i}$$ (11) 此外, 黏性应力张量$ \sigma_{ij} $, 则可以由非平衡态分布函数$ f_\alpha^{neq}\equiv f_\alpha-f_\alpha^{eq} $以及离散速度$ {\boldsymbol{\xi}}_\alpha $表示如下[11]
$$ \sigma_{ij}=-\sum\limits_{\alpha=1}^{q}f_\alpha^{neq} \xi_{\alpha,i}\xi_{\alpha,j} $$ (12) 而根据式 (11)和式(12) 以及Hermite多项式的具体表达式可知, 密度$ \rho $、速度$ {\boldsymbol{u}} $以及黏性应力张量$ {\boldsymbol{\sigma}} $正好分别对应连续分布函数变量$ f_\alpha $的第$ 0 $阶、第$ 1 $阶Hermite矩以及非平衡态分布函数$ f_\alpha^{neq} $的第$ 2 $阶Hermite矩
$$ \rho=a^{(0)},\qquad \rho u_i=a^{(1)}_i,\qquad \sigma_{ij}=-a^{neq,(2)}_{ij} $$ (13) 式中, ${\boldsymbol{a}}^{(n)}=\displaystyle\sum _{\alpha=1}^{q}f_\alpha {\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha)$; ${\boldsymbol{a}}^{neq,(n)}=\displaystyle\sum _{\alpha=1}^{q}f^{neq}_\alpha {\boldsymbol{H}}^{(n)} ({\boldsymbol{\xi}}_\alpha)$ 分别代表$ f_\alpha $和$ f_\alpha^{neq} $对应的第$ n $阶Hermite矩.
将式 (7) 代入到式 (5) 后, 可以反解出
$$ f_\alpha=g_\alpha-\frac{1}{2\tau_\nu}\left(g_\alpha -g_\alpha^{eq}+ \frac{\Delta t}{2}S_\alpha\right)+ \frac{\Delta t}{2}S_\alpha $$ (14) 式中, $ \tau_\nu $代表格子松弛系数, 其与网格分辨率有关
$$ \tau_\nu=\frac{\tau}{\Delta t }+ \frac{1}{2} $$ (15) 由于时空离散前后平衡态部分并没有发生任何改变, 故为了统一符号, 取$ g_\alpha^{eq}\equiv f_\alpha^{eq} $. 将式 (14) 和 式(15) 代入式 (7), 可以得到BGK碰撞模型下完全由新变量$ g_\alpha $表示的碰撞项$ \varOmega_{\alpha} $
$$ \varOmega_{{\rm{BGK}},\alpha}=-\frac{1}{\tau_\nu}\left(g_\alpha-g^{eq}_\alpha+\frac{\Delta t}{2}S_{\alpha}\right) $$ (16) 将式 (16) 代入式 (6), 便可以得到使用BGK碰撞模型的标准LBM演化方程
$$ \begin{split} &g_{\alpha}\left({\boldsymbol{x}}+{\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t\right)=\left(1-\frac{1}{\tau_\nu}\right)g_{\alpha}({\boldsymbol{x}}, t)+\\ &\qquad \frac{1}{\tau_\nu}g_\alpha^{eq}({\boldsymbol{x}}, t)+\Delta t\left(1-\frac{1}{2\tau_\nu}\right) S_{\alpha}({\boldsymbol{x}}, t)\end{split} $$ (17) 上述演化方程可以拆分为局部的碰撞和迁移两个过程. (1) 碰撞过程: $g_{\alpha}^\prime\left({\boldsymbol{x}}, t\right)=\left(1-\dfrac{1}{\tau_\nu}\right)g_{\alpha}({\boldsymbol{x}}, t)+ \dfrac{1}{\tau_\nu}g_\alpha^{eq}({\boldsymbol{x}}, t)+ \Delta t\left(1-\dfrac{1}{2\tau_\nu}\right) S_{\alpha}({\boldsymbol{x}}, t)$; (2) 迁移过程: $g_{\alpha}\left({\boldsymbol{x}}+ {\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t\right)=g^\prime_{\alpha}\left({\boldsymbol{x}}, t\right)$.
当$ S_\alpha $代表体力项$ F_\alpha=-{\boldsymbol{b}} \cdot \nabla_{{\boldsymbol{\xi}}_\alpha} f_\alpha $时, 其前两阶Hermite展开可以表示如下[28]
$$ \begin{split} &F_\alpha=w_\alpha \sum\limits_{n=0}^{2} \frac{1}{n!}{\boldsymbol{a}}^{(n)}_F({\boldsymbol{x}},t){\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha)=\\ &\qquad w_\alpha\rho\left[\underbrace{{\boldsymbol{b}}\cdot{\boldsymbol{\xi}}_\alpha}_{1{\rm{st}}}+\underbrace{({\boldsymbol{b}}\cdot{\boldsymbol{\xi}}_\alpha)({\boldsymbol{u}}\cdot{\boldsymbol{\xi}}_\alpha)-{\boldsymbol{b}}\cdot{\boldsymbol{u}}}_{2{\rm{nd}}\quad {\rm{order}}}\right]\end{split} $$ (18) 式中, $ {\boldsymbol{b}} $代表单位质量所受到的体力; ${\boldsymbol{a}}^{(n)}_F= \displaystyle\sum\limits_{\alpha=1}^{q}F_\alpha {\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha)$代表力项分布函数$ F_\alpha $对应的第$ n $阶Hermite矩, 其前两阶Hermite矩可以表示为
$$ a_F^{(0)}=0,\quad a^{(1)}_{F,i}=\rho b_i,\quad a^{(2)}_{F,ij}=\rho (u_i b_j+u_j b_i)$$ (19) 而通过Hermite展开得到的力项表达式 (18), 也与Guo等[32]通过反设计得到的结果完全一致.
除此之外, 由于式 (17) 的解$ g_\alpha\neq f_\alpha $, 故由离散分布函数变量$ g_\alpha $计算宏观量的表达式需要重新给出. 通过式 (14), 可以得到连续分布函数变量$ f_\alpha $对应的各阶Hermite矩$ {\boldsymbol{a}}^{(n)} $与离散分布函数变量$ g_\alpha $的各阶Hermite矩$ {\boldsymbol{a}}^{(n)}_g $间的关系
$$ {\boldsymbol{a}}^{(n)}={\boldsymbol{a}}^{(n)}_g-\frac{1}{2\tau_\nu}{\boldsymbol{a}}^{neq,(n)}_g+\frac{\Delta t}{2}\left(1-\frac{1}{2\tau_\nu}\right){\boldsymbol{a}}^{(n)}_s $$ (20) 其中
$$ \left.\begin{split} &{\boldsymbol{a}}^{(n)}_g=\sum\limits_{\alpha=1}^{q}g_\alpha {\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha),\quad {\boldsymbol{a}}^{neq,(n)}_g=\sum\limits_{\alpha=1}^{q}g^{neq}_\alpha {\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha)\\ &{\boldsymbol{a}}^{(n)}_s=\sum\limits_{\alpha=1}^{q}S_\alpha {\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha)\end{split}\right\} $$ (21) 且$ g_\alpha^{neq}\equiv g_\alpha-g_\alpha^{eq} $. 考虑$ f_\alpha^{eq}= g_\alpha^{eq} $后, 式 (20) 两边同时减去平衡态分布函数对应的各阶Hermite矩, 可以得到连续非平衡态分布函数变量$ f_\alpha^{neq} $的各阶Hermite矩$ {\boldsymbol{a}}^{neq,(n)} $与离散非平衡态分布函数变量$ g_\alpha^{neq} $的各阶Hermite矩$ {\boldsymbol{a}}_g^{neq,(n)} $间的关系如下
$$ {\boldsymbol{a}}^{neq,(n)}=\left(1-\frac{1}{2\tau_\nu}\right)\left({\boldsymbol{a}}_g^{neq,(n)}+\frac{\Delta t}{2}{\boldsymbol{a}}_s^{(n)}\right) $$ (22) 根据式 (22) 和式 (13), 可以显式表达出离散非平衡态分布函数变量$ g_\alpha^{neq} $的前两阶Hermite矩
$$ \left. \begin{split} &a^{neq,(0)}_g=-\frac{\Delta t}{2}a_s^{(0)}\\ &a^{neq,(1)}_{g,i}=-\frac{\Delta t}{2}a_{s,i}^{(1)}\\ &a^{neq,(2)}_{g,ij}=-\frac{2\tau_\nu}{2\tau_\nu-1}\sigma_{ij}-\frac{\Delta t}{2}a_{s,ij}^{(2)} \end{split} \right\} $$ (23) 当考虑源项为体力项时, 根据式 (19), $ {\boldsymbol{a}}_g^{neq,(n)} $的前两阶可以表示为
$$ \left. \begin{split} &a^{neq,(0)}_g=0\\ &a^{neq,(1)}_{g,i}=-\frac{\Delta t}{2}\rho b_i\\ &a^{neq,(2)}_{g,ij}=-\frac{2\tau_\nu}{2\tau_\nu-1}\sigma_{ij}-\frac{\Delta t}{2}\rho\left(b_i u_j+ b_j u_i\right) \end{split} \right\} $$ (24) 进而可以通过式 (13)、式(20) 以及式 (24) 得到当源项为体力项时, 由离散分布函数变量$ g_\alpha $表示的密度$ \rho $、速度$ {\boldsymbol{u}} $以及由离散非平衡态分布函数变量$ g_\alpha^{neq} $表示的黏性应力张量$ {\boldsymbol{\sigma}} $
$$ \rho=\sum\limits_{\alpha=1}^{q} g_\alpha,\quad \rho u_i=\sum\limits_{\alpha=1}^{q} \xi_{\alpha,i} g_\alpha+\frac{\Delta t}{2}\rho b_i \tag{25a}$$ $$ \sigma_{ij}=-\left(1-\frac{1}{2\tau_{\nu}}\right)\left[\sum\limits_{\alpha=1}^{q} \xi_{\alpha ,i}\xi_{\alpha ,j}g_{\alpha}^{neq}+\frac{\Delta t}{2}\rho(b_i u_j+ b_j u_i)\right] \tag{25b}$$ 最后, 在经过合适的初始化后, 式 (17) 便可以显式推进. 而当初始流场中包含非平衡态时, 正确的初始化分布函数对于一些非定常问题是至关重要的. 其中平衡态分布函数$ g_\alpha^{eq} $可以直接根据式 (8) 通过初始场的密度以及速度给出, 而非平衡态分布函数的初始化可以借由Hermite展开给出[33]. 对于等温的流动问题, 分布函数的非平衡态部分可以直接使用其前两阶Hermite展开来近似, 即
$$ g^{neq}_\alpha({\boldsymbol{x}},t)=w_\alpha\sum\limits_{n=0}^{2}\frac{1}{n!}{\boldsymbol{a}}_g^{neq,(n)}({\boldsymbol{x}},t){\boldsymbol{H}}^{(n)}({\boldsymbol{\xi}}_\alpha) $$ (26) 将式 (24) 代入到上式中, 便可以得到初始流场在BGK模型下且考虑源项为体力项时, 非平衡态分布函数$ g_\alpha^{neq} $的具体表达形式
$$ g_\alpha^{neq}=w_\alpha \left\{\underbrace{-\frac{\Delta t}{2}\rho b_i\xi_{\alpha,i}}_{1{\rm{st}}}\underbrace{-\frac{1}{2}\left[\frac{2\tau_\nu}{2\tau_\nu-1}(\sigma_{ij}\xi_{\alpha,i}\xi_{\alpha,j}-\sigma_{ii})+\Delta t \rho (b_i\xi_{\alpha,i}u_j\xi_{\alpha,j}-b_iu_i)\right]}_{2{\rm{nd}}\ {\rm{order}}}\right\} $$ (27) 由于我们的目的是保证初始时刻分布函数所对应的密度、速度以及黏性应力张量与初始条件一致, 故这里只是展到二阶.
以上便是从DVBE到LBM的详细推导过程. 其中注意的是, 在离散过程中引入了离散分布函数变量 $ g_\alpha $, 而且核心演化关系式 (17) 实际求解的也是离散分布函数变量$ g_\alpha $, 而不是时空离散前的连续分布函数变量$ f_\alpha $. 此外, 从式 (5) 不难看出, 离散分布函数变量$ g_\alpha $与时间步长$ \Delta t $相关, 进而也与网格分辨率有关. 故当使用式 (17) 作为计算流域上的演化公式时, 即使忽略时空离散误差, 同一时空点下不同网格分辨率对应的$ g_\alpha $仍然是不一致的, 所以粗细网格重合点上的$ g_\alpha $并不能直接交换, 而这也是局部网格加密算法中不同网格分辨率在重合点处离散分布函数$ g_\alpha $需要转换的原因[16].
1.2 MRT碰撞模型
MRT模型是d'Humieres [24]在1992年提出的一种碰撞模型, 该模型可以被视为广义的BGK模型. 相比于只有一个松弛参数的原始BGK模型, MRT模型中则有着多个松弛参数用以描述不同物理量趋向平衡态的过程. 除了少数几个有着明确物理含义的松弛参数外, MRT模型还有若干个可调参数, 相关文献[34-35]在对MRT模型进行系统研究后指出, 对这些参数进行合理的调节可以有效改善数值稳定性以及减小边界条件的误差.
MRT模型对应的速度离散后的碰撞项可以表示如下[36]
$$ |\varOmega_{{\rm{MRT}}}({\boldsymbol{x}},t)\rangle=-{\boldsymbol{\varLambda}}\left(|f({\boldsymbol{x}},t)\rangle-|f^{eq} ({\boldsymbol{x}},t)\rangle\right) $$ (28) 其中$ {\boldsymbol{\varLambda}} $代表物理松弛系数矩阵, 其与网格分辨率无关; 而$ f^{eq}_\alpha $的表达式仍然使用式 (8). 上式中使用了与D'Humieres相同的记号, 以“$ \langle\cdot| $”, “$ |\cdot\rangle $”来分别表示行向量以及列向量, 即
$$ \begin{split} &|f({\boldsymbol{x}},t)\rangle\equiv\left(f_1,f_2,\cdots,f_{q}\right)^{\rm{T}},\; \; |f^{eq} ({\boldsymbol{x}},t)\rangle\equiv\left(f^{eq}_1,f^{eq}_2,\cdots,f^{eq}_{q}\right)^{\rm{T}}\\ &|\varOmega_{{\rm{MRT}}}({\boldsymbol{x}},t)\rangle \equiv\left(\varOmega_{{\rm{MRT}},1},{{\varOmega}}_{{\rm{MRT}},2},\cdots,\varOmega_{{\rm{MRT}},{q}}\right)^{\rm{T}}\end{split}$$ 这里上标“${\rm{T}}$”代表转置操作符. 将式 (28) 代入到式 (5) 的矩阵形式中后, 可以反解出由离散分布函数变量$ |g\rangle $所表示的连续分布函数变量$ |{{f}}\rangle $如下
$$ |f({\boldsymbol{x}},t)\rangle=\left({\boldsymbol{I}}+\frac{\Delta t}{2}{\boldsymbol{\varLambda}}\right)^{-1}\left(|g\rangle-|g^{eq}\rangle+\frac{\Delta t}{2}|S\rangle\right) + |g^{eq}\rangle $$ (29) 其中$ {\boldsymbol{I}} $为单位矩阵; 上标“$ -1 $”代表矩阵求逆操作符. 进而可以得到MRT模型下由离散分布函数变量$ |g\rangle $表示的碰撞项
$$ |\varOmega_{{\rm{MRT}}}({\boldsymbol{x}},t)\rangle=-{\boldsymbol{\varLambda}}\left({\boldsymbol{I}}+\frac{\Delta t}{2}{\boldsymbol{\varLambda}}\right)^{-1}\left(|g\rangle-|g^{eq}\rangle+\frac{\Delta t}{2}|S\rangle\right) $$ (30) 将式 (30) 代入到式 (6) 的矩阵形式中, 可以得到
$$ \begin{split} &|g\left({\boldsymbol{x}}+{\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t\right)\rangle=|g({\boldsymbol{x}}, t)\rangle-\\ &\qquad \Delta t {\boldsymbol{\varLambda}}\left({\boldsymbol{I}}+\frac{\Delta t}{2}{\boldsymbol{\varLambda}}\right)^{-1}\left(|g({\boldsymbol{x}}, t)\rangle-|g^{eq}({\boldsymbol{x}}, t)\rangle\right)+\\ &\qquad \Delta t \left[{\boldsymbol{I}}-\frac{1}{2}\Delta t {\boldsymbol{\varLambda}}\left({\boldsymbol{I}}+\frac{\Delta t}{2}{\boldsymbol{\varLambda}}\right)^{-1}\right]|S({\boldsymbol{x}}, t)\rangle \end{split} $$ (31) 若记
$$ \Delta t{\boldsymbol{\varLambda}}\left({\boldsymbol{I}}+\frac{\Delta t}{2}{\boldsymbol{\varLambda}}\right)^{-1}={\boldsymbol{M}}^{-1}{\boldsymbol{\hat{S}}}{\boldsymbol{M}} $$ (32) 则将上式代入到式 (31) 中后, 便可以得到MRT模型下关于离散分布函数变量$ |g\rangle $的演化方程
$$ \begin{split} &|g\left({\boldsymbol{x}}+{\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t\right)\rangle=|g({\boldsymbol{x}}, t)\rangle-\\ &\qquad {\boldsymbol{M}}^{-1}{\boldsymbol{\hat{S}}}{\boldsymbol{M}}\left(|g({\boldsymbol{x}}, t)\rangle-|g^{eq}({\boldsymbol{x}}, t)\rangle\right)+\\ &\qquad \Delta t \left({\boldsymbol{I}}-\frac{1}{2}{\boldsymbol{M}}^{-1}{\boldsymbol{\hat{S}}}{\boldsymbol{M}}\right)|S({\boldsymbol{x}}, t)\rangle \end{split} $$ (33) 其中, $ {\boldsymbol{M}} $是一个$ q\times q $的正交变换矩阵; $ {\boldsymbol{M}}^{-1} $为$ {\boldsymbol{M}} $的逆矩阵, 且满足以下关系
$$ |m({\boldsymbol{x}}, t)\rangle={\boldsymbol{M}}|g({\boldsymbol{x}}, t)\rangle,\; |g({\boldsymbol{x}}, t)\rangle={\boldsymbol{M}}^{-1}|m({\boldsymbol{x}}, t)\rangle $$ (34) 式中, $ |m\rangle $代表转换矩阵$ {\boldsymbol{M}} $下分布函数$ |g\rangle $对应的一组矩. 对式 (33) 两边同时左乘转换矩阵$ {\boldsymbol{M}} $后便可得到MRT模型下关于$ |m\rangle $的演化方程[23]
$$ \begin{split} &|m({\boldsymbol{x}}+{\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t)\rangle=|m({\boldsymbol{x}}, t)\rangle-\\ &\qquad{\boldsymbol{\hat{S}}}\left(|m({\boldsymbol{x}}, t)\rangle- |m^{eq}({\boldsymbol{x}}, t)\rangle\right)+\Delta t\left({\boldsymbol{I}}-\frac{1}{2}{\boldsymbol{\hat{S}}}\right)|m_{s}({\boldsymbol{x}}, t)\rangle\end{split} $$ (35) 这里$ |m_s({\boldsymbol{x}},t)\rangle={\boldsymbol{M}}|S({\boldsymbol{x}},t)\rangle $. 而转换矩阵$ {\boldsymbol{M}} $与速度集的选取有关, 以D2Q9速度集为例, 其对应的转换矩阵$ {\boldsymbol{M}} $为 [34]
$$ {\boldsymbol{M}}=\left[\begin{array}{ccccccccc} 1&1&1&1&1&1&1&1&1 \\ -4&-1&-1&-1&-1&2&2&2&2 \\ 4&-2&-2&-2&-2&1&1&1&1 \\ 0&1&0&-1&0&1&-1&-1&1 \\ 0&-2&0&2&0&1&-1&-1&1 \\ 0&0&1&0&-1&1&1&-1&-1 \\ 0&0&-2&0&2&1&1&-1&-1 \\ 0&1&-1&1&-1&0&0&0&0 \\ 0&0&0&0&0&1&-1&1&-1 \end{array}\right] $$ (36) 而格子松弛系数矩阵$ {\boldsymbol{\hat{S}}} $为一个对角矩阵, 其对角线上包含着各个矩演化过程的松弛参数, 且由式 (32) 可知其与网格分辨率有关. 对于D2Q9速度集而言, 该矩阵可以表示如下
$$ {\boldsymbol{\hat{S}}}={\rm{diag}}\left[0, \ s_e, \ s_\epsilon, \ 0, \ s_q, \ 0, \ s_q, \ s_\nu, \ s_\nu\right] $$ (37) 其中, $ s_\nu $与运动黏度系数$ \nu $, 以及$ s_e $与体积黏度系数$ \zeta $的关系可以由CE分析给出[37]
$$ \nu=\left(\frac{1}{s_\nu}-\frac{1}{2}\right)\Delta t,\quad \zeta =\left(\frac{1}{s_e}-\frac{1}{2}\right)\Delta t $$ (38) 而BGK模型则相当于${\boldsymbol{\hat{S}}}_{{\rm{BGK}}}=\dfrac{1}{\tau_\nu}{\boldsymbol{I}}$, 此时MRT模型下分布函数的演化方程式 (33) 将退化为BGK模型下分布函数的演化方程式 (17).
与BGK模型类似, MRT模型下演化关系式 (33) 的解$ |g\rangle\neq |f\rangle $, 故也需要给出MRT模型下, 由离散分布函数变量$ |g\rangle $计算宏观量的表达式. 通过引入Hermite多项式矩阵$ {\boldsymbol{H}} $和Hermite转换矩阵$ {\boldsymbol{M}}_H $, 可以根据式 (29) 得到MRT模型下, 连续分布函数变量$ |f\rangle $对应的各阶Hermite矩$ |a\rangle $与离散分布函数变量$ |g\rangle $对应的各阶Hermite矩$ |a_g\rangle $间的关系; 以及连续非平衡态分布函数变量$ |f^{neq}\rangle $对应的各阶Hermite矩$ |a^{neq}\rangle $与离散非平衡态分布函数变量$ |g^{neq}\rangle $对应的各阶Hermite矩$ |a_g^{neq}\rangle $间的关系如下(各变量的具体表达式以及详细推导过程见附录A)
$$ |{{a}}\rangle={\boldsymbol{M}}_H^{-1}\left(|a_g\rangle+\frac{\Delta t}{2} |a_s\rangle\right)+\left({\boldsymbol{I}}-{\boldsymbol{M}}_H^{-1}\right)|a_g^{eq}\rangle \tag{39a}$$ $$ |{{a}}^{neq}\rangle={\boldsymbol{M}}_H^{-1}\left(|a_g^{neq}\rangle+\frac{\Delta t}{2} |a_s\rangle\right) \tag{39b}$$ 进而我们可以通过式 (13)、式(19)、式(39 a)以及附录A中MH得到MRT模型下当源项为体力项时, 由离散分布函数变量$ g_\alpha $表示的密度$ \rho $、速度$ {\boldsymbol{u}} $
$$ \rho=\sum\limits_{\alpha=1}^{q} g_\alpha,\quad \rho u_i=\sum\limits_{\alpha=1}^{q} \xi_{\alpha,i} g_\alpha+\frac{\Delta t}{2}\rho b_i $$ (40) 以及通过附录A中MH中参数表达式和$a_g^{neq} $的各阶Hermite矩显示, 可得到MRT模型下, 由离散非平衡态分布函数变量$ g_\alpha^{neq} $的Hermite矩和力项分布函数$ F_\alpha $的Hermite矩表示的黏性应力张量$ {\boldsymbol{\sigma}} $
$$ \begin{split} &\sigma_{x x} = -\left(a_{g, x x}^{{neq, (2) }} + \frac{\Delta t}{2} a_{F, x x}^{(2)}\right)+\frac{1}{4} s_e\biggl[\left(a_{g, x x}^{n e q,(2)}+a_{g, y y}^{n e q,(2)}\right)+\\ &\qquad \frac{\Delta t}{2}\left(a_{F, x x}^{(2)}+a_{F, y y}^{(2)}\right)\biggl] +\frac{1}{4} s_v\biggl[\left(a_{g, x x}^{n e q,(2)}-a_{g, y y}^{n e q,(2)}\right)+\\ &\qquad \frac{\Delta t}{2}\left(a_{F, x x}^{(2)}-a_{F, y y}^{(2)}\right)\biggl]\end{split}\tag{41a} $$ $$ \begin{split} &\sigma_{y y}-\left(a_{g, y y}^{n e q,(2)}+\frac{\Delta t}{2} a_{F, y y}^{(2)}\right)+\frac{1}{4} s_e\biggl[\left(a_{g, x x}^{n e q,(2)}+a_{g, y y}^{n e q,(2)}\right)+\\ &\qquad \frac{\Delta t}{2}\left(a_{F, x x}^{(2)}+a_{F, y y}^{(2)}\right)\biggl]-\frac{1}{4} s_v\biggl[\left(a_{g, x x}^{n e q,(2)}-a_{g, y y}^{n e q,(2)}\right)+\\ &\qquad \frac{\Delta t}{2}\left(a_{F, x x}^{(2)}-a_{F, y y}^{(2)}\right)\biggl]\end{split}\tag{41b} $$ $$ \sigma_{xy}=-\left(1-\frac{1}{2}s_\nu\right)\left(a_{g, x y}^{n e q,(2)}+\frac{\Delta t}{2}a_{F, x y}^{(2)}\right) \tag{41c}$$ 上述结果也与Chen等[37]通过CE分析得到结果保持一致.
当初始流场中包含非平衡态时, MRT模型下的分布函数也需要正确的初始化. 同样可以使用Hermite 展开来初始化MRT模型下的分布函数, 其中平衡态分布函数的初始化与BGK模型一致, 而离散非平衡态分布函数$ g_\alpha^{neq} $的初始化, 则需要将附录A中$|g^{neq} \rangle$的前6个矩代入到式 (26), 进而可以得到MRT模型下当源项为体力项时, 初始流场中$ g_\alpha^{neq} $的具体表达式
$$ g_\alpha^{neq}({\boldsymbol{x}},t)=w_\alpha \left[ \begin{split} &\underbrace{-\frac{\Delta t}{2}\rho\left(b_{x}\xi_{\alpha,x}+b_{y}\xi_{\alpha,y}\right)}_{1{\rm{st}}\ {\rm{order}}}\underbrace{-\frac{1}{2C_2}(c_{22}\sigma_{xx}+c_{23}\sigma_{yy})(\xi_{\alpha ,x}\xi_{\alpha ,x}-1)}_{2{\rm{nd}}\ {\rm{order}}}\\ &\underbrace{-\frac{1}{2}d_2\sigma_{xy}\xi_{\alpha ,x}\xi_{\alpha ,y}-\frac{1}{2C_2}(c_{23}\sigma_{xx}+c_{22}\sigma_{yy})(\xi_{\alpha ,y}\xi_{\alpha ,y}-1)}_{2{\rm{nd}}\ {\rm{order}}}\\ &\underbrace{-\frac{\Delta t}{2} \rho (b_i\xi_{\alpha,i}u_j\xi_{\alpha,j}-b_iu_i)}_{2{\rm{nd}}\ {\rm{order}}} \end{split} \right] $$ (42) 上式中各参数的具体形式可以参见附录A.
此外, 由于矩阵$ {\boldsymbol{M}} $和矩阵$ {\boldsymbol{\hat{S}}} $初始化后就已确定, 且其在式 (33) 演化过程中不再改变, 故可以在初始化时将$ {\boldsymbol{M}}^{-1}{\boldsymbol{\hat{S}}}{\boldsymbol{M}} $计算后单独存储起来以避免在执行碰撞步骤时重复计算, 从而节省计算量. 且直接使用式 (33) 去执行分布函数的碰撞步骤相比于使用式 (34) 和式 (35) 更加节省时间, 经我们程序测试后发现CPU消耗时间减少$25\text{%}$左右.
2. 局部网格加密中的数据转换
如前所述, 同一点处不同网格分辨率下的离散分布函数变量$ g_\alpha $之间的转换是必要且不能忽略的. 使用局部网格加密算法需要保证, 当粗细网格均使用同一离散速度集且在忽略时空离散误差后, 粗细网格区域在同一时空点下数值求解后得到的连续分布函数变量$ f_\alpha $必须是一致的, 而且不同分辨率下同一点处的物性参数如: 物理松弛系数$ \tau $以及运动黏度系数$ \nu $等也必须保持一致.
LBM在对式 (1) 做时空离散时为了消除隐式, 引入了与分辨率有关的离散分布函数变量$ g_\alpha $以及相应的格子松弛系数$ \tau_\nu $(BGK模型)和$ \hat{{\boldsymbol{S}}} $(MRT模型)
$$ \left.\begin{split} & g_\alpha({\boldsymbol{x}},t)=f_\alpha({\boldsymbol{x}},t) - \frac{\Delta t}{2}(\varOmega_{\alpha}({\boldsymbol{x}},t) + S_{\alpha}({\boldsymbol{x}},t))\\ & \tau_\nu=\frac{\tau}{\Delta t}+\frac{1}{2} \quad \text{(BGK model)}\\ &|g({\boldsymbol{x}},t)\rangle=|f({\boldsymbol{x}},t)\rangle - \frac{\Delta t}{2}\left(|\varOmega({\boldsymbol{x}},t)\rangle + |S({\boldsymbol{x}},t)\rangle\right)\\ & \hat{{\boldsymbol{S}}}=\Delta t{\boldsymbol{M}} {\boldsymbol{\varLambda}}\left({\boldsymbol{I}}+\frac{\Delta t}{2}{\boldsymbol{\varLambda}}\right)^{-1}{\boldsymbol{M}}^{-1}\; \quad \text{(MRT model)} \end{split}\right\} $$ (43) 根据上式可知, 不同网格分辨率下的同一时空点上, 格子松弛系数$ \tau_\nu $和$ \hat{\boldsymbol {S}} $以及实际求解得到的离散分布函数变量$ g_\alpha $和$ |g\rangle $是不一致的, 故粗细网格在交换这些量之前需要经过合适的转换. 与之相反的是, 物理松弛系数$ \tau $以及$ {\boldsymbol{\varLambda}}^{-1} $是不随网格分辨率变化的; 而当考虑粗细网格区域均使用同一离散速度集时, 则数值求解得到的连续分布函数变量$ f_\alpha $和$ |f\rangle $, 在不考虑时空离散误差的前提下, 也是不随网格分辨率变化的. 故接下来将以此为桥梁, 建立起BGK模型和MRT模型下, 不同网格分辨率间的数据转换关系.
下文中, 与粗网格和细网格区域相关的任何量分别用上标“$ c $”和“$ f $”表示. 设局部加密时粗细网格的空间步长比例为
$$ n=\Delta x^c/\Delta x^f $$ (44) 局部网格加密算法中不同网格分辨率下的区域仍然使用的是同一离散速度集, 故在考虑对流尺度下, 粗细网格的时间步长也有同样的比例关系
$$ n=\Delta t^c/\Delta t^f $$ (45) 故为了保证粗细网格上物理时间的统一推进, 粗网格上运行一步则细网格上需要运行$ n $步. 故$ n $常设为大于$ 1 $的整数.
2.1 BGK 模型下的局部网格加密
为了将我们的方法与现有推导过程更好地对比, 在给出新的方法前, 先梳理一下Liou等[25]在考虑源项时, 分布函数转换关系的推导过程.
首先对式 (17) 使用泰勒展开
$$ \begin{split} &\Delta t\left(\frac{\partial}{\partial t}+{\boldsymbol{\xi}}_\alpha\cdot\nabla\right)g_\alpha+\frac{\Delta t^2}{2}\left(\frac{\partial}{\partial t}+{\boldsymbol{\xi}}_\alpha\cdot\nabla\right)^2g_\alpha=\\ &\qquad -\frac{1}{\tau_\nu}\left(g_{\alpha}-g_{\alpha}^{eq}\right)+\left(1-\frac{1}{2\tau_\nu}\right)\Delta t S_{\alpha}\end{split} $$ (46) 考虑CE展开如下
$$ g_\alpha=g_\alpha^{(0)} + \epsilon g_\alpha^{(1)}+ \epsilon^2 g_\alpha^{(2)}+\cdots \tag{47a} $$ $$ \frac{\partial}{\partial _t}=\epsilon\frac{\partial}{\partial{t_1}}+\epsilon^2\frac{\partial}{\partial{t_2}}\tag{47b}$$ $$ \nabla=\epsilon\nabla_1,\ S_\alpha=\epsilon S^{(1)}_\alpha \tag{47c} $$ 其中$ \epsilon $为与克努森数$ Kn $同阶的小量. 将式 (47) 代入到式 (46) 中, 便可以得到与$ \epsilon^0 $,$ \epsilon^1 $,$ \epsilon^2 $同阶的方程
$$ \epsilon^0: g_\alpha^{(0)}=g_\alpha^{eq}\tag{48a} $$ $$ \epsilon^1:\left(\frac{\partial}{\partial{t_1}}+{\boldsymbol{\xi}}_\alpha\cdot\nabla_1\right)g_\alpha^{(0)}-S_\alpha^{(1)}=-\frac{1}{\tau_\nu\Delta t}\left(g_\alpha^{(1)}+\frac{\Delta t}{2}S_\alpha^{(1)}\right)\tag{48b} $$ $$ \begin{split} &\epsilon^2: \frac{\partial g_\alpha^{(0)}}{\partial t_2}+\left(1-\frac{1}{2\tau_\nu}\right)\left(\frac{\partial}{\partial{t_1}}+{\boldsymbol{\xi}}_\alpha\cdot\nabla_1\right)\cdot\\ &\qquad \left(g_\alpha^{(1)}+\frac{\Delta t}{2}S_\alpha^{(1)}\right)=-\frac{1}{\tau_\nu\Delta t}g_\alpha^{(2)}\end{split}\tag{48c} $$ 而式 (48b) 等号左边项均为宏观量与离散速度$ {\boldsymbol{\xi}}_\alpha $的函数, 故当粗细网格区域使用同一离散速度集时, 其与网格分辨率无关, 故如下关系式成立
$$ \frac{1}{\tau_\nu^c\Delta t^c}\left(g_\alpha^{(1),c}+\frac{\Delta t^c}{2}S_\alpha^{(1)}\right)=\frac{1}{\tau_\nu^f\Delta t^f}\left(g_\alpha^{(1),f}+\frac{\Delta t^f}{2}S_\alpha^{(1)}\right) $$ (49) 使用$ \epsilon g_\alpha^{(1)} \approx g_\alpha^{neq} $假设后, 则可以得到离散非平衡态分布函数变量$ g_\alpha^{neq} $的转换关系如下
$$ \frac{1}{\tau_\nu^c\Delta t^c}\left(g_\alpha^{neq,c}+\frac{\Delta t^c}{2}S_\alpha\right)=\frac{1}{\tau_\nu^f\Delta t^f}\left(g_\alpha^{neq,f}+\frac{\Delta t^f}{2}S_\alpha\right) $$ (50) 而格子松弛系数$ \tau_\nu $的转换则是通过保证动力黏度系数 $ \mu=p(\tau_\nu-0.5)\Delta t $ 一致的基础上推导出来的, 即
$$ p\left(\tau_\nu^c-\frac{1}{2}\right)\Delta t^c=p\left(\tau_\nu^f-\frac{1}{2}\right)\Delta t^f $$ (51) 以上便是Liou等在考虑源项存在时, 给出的BGK模型下粗细网格间离散非平衡态分布函数变量$ g_\alpha^{neq} $以及格子松弛系数$ \tau_\nu $转换关系的推导过程, 而$ g_\alpha^{neq} $转换关系的推导过程本质上是利用CE展开寻找与网格分辨率无关的量. 不难发现整个推导过程是较为复杂的, 而且使用了$ \epsilon g_\alpha^{(1)} $来近似$ g_\alpha^{neq} $, 对于恢复N-S方程而言, 这样的近似是足够的, 但对于需要捕捉稀薄效应的高阶LBM方法 [38]而言, 一阶近似显然是不够的. 下面将给出一种更加简洁且直观的推导方法.
如前所述, BGK模型下的物理松弛系数$ \tau $ 和 连续分布函数变量$ f_\alpha $与网格分辨率无关, 故由式 (15) 和式 (14) 可知如下关系是成立的
$$ \left(\tau^c_\nu-\frac{1}{2}\right)\Delta t^c=\left(\tau^f_\nu-\frac{1}{2}\right)\Delta t^f\tag{52a} $$ $$\begin{split} &g^c_\alpha+\frac{1}{2\tau^c_\nu}(g_\alpha^{eq}-g^c_\alpha)+\frac{(2\tau^c_\nu-1)\Delta t^c}{4\tau^c_\nu}S_{\alpha}=\\ &\qquad g^f_\alpha+\frac{1}{2\tau^f_\nu}(g_\alpha^{eq}-g^f_\alpha)+\frac{(2\tau^f_\nu-1)\Delta t^f}{4\tau^f_\nu}S_{\alpha}\end{split}\tag{52b} $$ 对上式简单整理后, 便可以得到不同网格分辨率在重合点处的数据转换关系如下
$$ \tau_\nu^f =\frac{1}{2} + n(\tau_\nu^c -\frac{1}{2})\tag{53a} $$ $$ g_\alpha^c =g_\alpha^f+\frac{(n-1)}{2\tau_\nu^f}\left(g_\alpha^f -g_\alpha^{eq}+\frac{\Delta t^f}{2}S_{\alpha}\right)-\frac{(n-1)}{2}\Delta t^f S_{\alpha}\tag{53b} $$ $$ g_\alpha^f =g_\alpha^{c}-\frac{(n-1)}{2n \tau_{v}^{c}}\left(g_\alpha^{c}-g_\alpha^{e q}+\frac{\Delta t^c}{2}S_{\alpha}\right)+\frac{(n-1)}{2n}\Delta t^cS_{\alpha}\tag{53c} $$ 将式 (52b) 二边减去平衡态并利用式 (52a), 便可以得到与式 (50) 一致的转换关系式.
相比之下, 我们的推导过程中除忽略时空离散误差外, 并没有做任何额外的假设(没有对非平衡态分布函数做近似假设), 故上述转换关系同样保证了粗细网格重合点处高阶非平衡态分布函数的一致. 然而, Liou等在使用$ \epsilon g_\alpha^{(1)} \approx g_\alpha^{neq} $的假设后, 也得到了与我们相同的结果. 而这主要是由于各阶非平衡态分布函数间的递推关系导致的(详细推导过程见附录 B)
$$ f_\alpha^{neq}=f_\alpha-f_\alpha^{eq}=\sum\limits_{k=1}^{N}\tau^k f_\alpha^{(k)}+O(\tau^{N+1}) \tag{54a} $$ $$ f_\alpha^{(k)}=-\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha^{(k-1)}=(-1)^{k-1}\left(\frac{{\rm{D}}}{{\rm{D}}t}\right)^{k-1}f_\alpha^{(1)} \tag{54b}$$ 从上式不难发现由于递推关系的存在, 实际上只要保证了$ f_\alpha^{(1)} $及$ \tau $在不同网格分辨下的一致性, 也就保证了其余阶数的一致性, 进而保证了$ f_\alpha^{neq} $的一致性. 而Liou等也恰恰是通过式 (49) 保证了不同分辨率下$ f_\alpha^{(1)} $的一致, 而关于不同分辨率下$ \tau $的一致性则是通过保证$ \mu $的一致来达到. 这也就是为何其只确保了$ f_\alpha^{(1)} $的一致, 但却可以得到一个更为“普适”的一致性.
上述推导过程不仅解释了我们的推导结果与前人结果一致的原因, 同时也为之前学者对非平衡态分布函数做一阶近似提供了理论依据, 而且上述推导过程也表明了局部网格加密算法同样可以被扩展到高阶LBM中.
2.2 MRT模型下的局部网格加密
下面将继续使用新方法给出MRT模型下, 不同网格分辨率在重合点处离散分布函数变量$ |g^c\rangle $和$ |g^f\rangle $以及格子松弛系数矩阵$ {\boldsymbol{\hat{S}}}^c $和$ {\boldsymbol{\hat{S}}}^f $间的转换关系. 与上一节BGK模型下的推导过程类似, 不同分辨率下的物理松弛系数矩阵$ {\boldsymbol{\varLambda}}^{-1} $应保持一致, 故由附录A可得
$$ \Delta t^c({\boldsymbol{M}}^{-1} {\boldsymbol{\hat{S}}}^c{\boldsymbol{M}})^{-1}-\frac{\Delta t^c}{2}{\boldsymbol{I}}=\Delta t^f({\boldsymbol{M}}^{-1} {\boldsymbol{\hat{S}}}^f{\boldsymbol{M}})^{-1}-\frac{\Delta t^f}{2}{\boldsymbol{I}} $$ (55) 对上式做简单整理后, 可得粗细网格间格子松弛系数矩阵$ \hat{{\boldsymbol{S}}} $中各元素间的关系
$$ \frac{1}{s_\beta^f}=n\left(\frac{1}{s_\beta^c}-\frac{1}{2}\right)+\frac{1}{2}\qquad(s_\beta^c\neq 0;s_\beta^c\neq 0;\beta=1,2,\cdots, q) $$ (56) 这里$ s_\beta^f $, $ s_\beta^c $分别代表对角阵$ {\boldsymbol{\hat{S}}}^f $和$ {\boldsymbol{\hat{S}}}^c $中的非零元素, 以D2Q9为例
$$ \frac{1}{s_\nu^f}=n\left(\frac{1}{s_\nu^c}-\frac{1}{2}\right)+\frac{1}{2},\qquad \frac{1}{s_e^f}=n\left(\frac{1}{s_e^c}-\frac{1}{2}\right)+\frac{1}{2} $$ (57) $$ \frac{1}{s_\epsilon^f}=n\left(\frac{1}{s_\epsilon^c}-\frac{1}{2}\right)+\frac{1}{2},\qquad \frac{1}{s_q^f}=n\left(\frac{1}{s_q^c}-\frac{1}{2}\right)+\frac{1}{2} $$ (58) 而由CE分析可知, $ s_q $和$ s_\epsilon $并不会直接出现在N-S方程中, 故可以被视为自由参数, 但其对数值稳定性及边界条件误差有一定影响 [34]. Chen等 [37]认为这两个参数在不同网格分辨率下可以不需要转换, 即$ s_\epsilon^f=s_\epsilon^c;\ s_q^f=s_q^c $. 但本文中为了保持转换关系的一致性, 仍然采用式 (58) 作为粗细网格间$ s_\epsilon $和$ s_q $的转换关系式.
同样, 直接从式 (29) 出发, 基于连续分布函数变量$ |f\rangle $在不同网格分辨率下的重合点处保持一致, 便可以得到MRT模型中粗细网格在时空重合点处的离散分布函数变量$ |g\rangle $的转换关系
$$\begin{split} &|g^{c}\rangle=|g^{f}\rangle+\frac{(n-1)}{2}({\boldsymbol{M}}^{-1} {\boldsymbol{\hat{S}}}_f{\boldsymbol{M}})\cdot\\ &\qquad \left(|g^{f}\rangle-|g^{eq}\rangle+\frac{\Delta t^f }{2}{\boldsymbol{S}}\right)-\frac{(n-1)}{2}\Delta t^f {\boldsymbol{S}}\end{split} \tag{59a}$$ $$ \begin{split} &|g^{f}\rangle=|g^{c}\rangle-\frac{(n-1)}{2n}({\boldsymbol{M}}^{-1} {\boldsymbol{\hat{S}}}_c{\boldsymbol{M}})\cdot\\ &\qquad\left(|g^{c}\rangle-|g^{eq}\rangle+\frac{\Delta t^c }{2}{\boldsymbol{S}}\right)+\frac{(n-1)}{2n}\Delta t^c {\boldsymbol{S}}\end{split}\tag{59b} $$ 且该结果与Huang等[23]的结果也完全一致.
综上, 我们以保证不同网格分辨间连续分布函数变量$ f_\alpha $, $ |f\rangle $以及物理松弛系数$ \tau $, $ {\boldsymbol{\varLambda}}^{-1} $一致为基础, 以更加直观且简洁的方法分别推导了BGK模型和MRT模型下, 不同网格分辨率在时空重合点处离散分布函数变量$ g_\alpha $, $ |g\rangle $和格子松弛系数$ \tau_\nu $, $ \hat{{\boldsymbol{S}}} $间的转换关系, 其相比于之前学者[23,25]的推导过程更加简洁且直观.
3. 局部网格加密算法的数值验证
本节数值模拟中将采用与Yu等[18]相同的网格排布方式, 其示意图如图2 所示. 同时, 为了避免局部网格加密算法中插值误差对数值模拟的影响, 我们使用与Gendre等[39]相同的高阶插值公式. 而对于加密界面处缺失的离散分布函数, 采用Chen等[40]提出的只使用空间插值的方法.
![]()
图 2 粗细网格排布示意图. 灰色区域为粗细网格的重叠区域Figure 2. Schematic diagram of the coarse and fine grid arrangement. The grey area is the overlap between the coarse and fine grids下面将通过使用局部网格加密技术的LBM来对强迫泰勒−格林涡流动、平板泊肃叶流中对流−扩散问题、顶盖驱动方腔流动和一维剪切波问题进行模拟, 在前两个算例中均包含非均匀且非定常的源项, 故以此来验证粗细网格间数据转换关系对复杂源项的适应性; 而有着较为复杂流动现象的顶盖驱动方腔流动则可以较好地验证局部网格加密技术的有效性; 一维剪切波问题的模拟则是用来检验局部网格加密技术对原本的LBM算法在耗散方面的影响.
3.1 强迫泰勒−格林涡流动
在标准的 Taylor-Green vortex 流动中没有外力的参与, 因此流动会随时间单调衰减. 为了验证源项存在时, 粗细网格间数据转换关系的正确性, 我们以Min等[41]使用的二维强迫泰勒−格林涡流动作为验证算例. 该流动中不仅有体力的参与, 且该体力非均匀且非定常, 故可以较为充分地验证BGK模型以及MRT模型下转换关系的正确性.
该流动问题对应的控制方程为
$$ \frac{\partial {\boldsymbol{u}}}{\partial t} + {\boldsymbol{u}}\cdot \nabla {\boldsymbol{u}}=-\frac{1}{\rho}\nabla p + \nu \nabla^2{\boldsymbol{u}}+{\boldsymbol{b}} $$ (60) 流动的物理区域被设为: $ (x,y)\in (0:L_x,0:L_y) $, 4周为周期边界条件, 考虑单位质量受到的体力为 $ {\boldsymbol{b}}=k^2\nu(1-Q){\boldsymbol{u}} $时, 则该流动对应的理论解如下
$$ \left.\begin{split} &u_x(x,y,t)=-U_0\cos(k_x x)\sin(k_y y){\rm{e}}^{-Qk^2\nu t}\\ &u_y(x,y,t)=\frac{k_x}{k_y}U_0\cos(k_y y)\sin(k_x x){\rm{e}}^{-Qk^2\nu t}\\ &p(x,y,t)=-\frac{1}{4}U_0^2\left[\cos(2k_x x)+\left(\frac{k_x}{k_y}\right)^2{\rm{cos}}(2k_y y)\right]\cdot \\ &\qquad {\rm{e}}^{-2Qk^2\nu t} + P_0\\ &\sigma_{xx}(x,y,t)=-\sigma_{yy}(x,y,t)=2\mu k_x U_0\sin(k_x x)\cdot \\ &\qquad \sin(k_y y){\rm{e}}^{-Qk^2\nu t}\\ &\sigma_{xy}(x,y,t)=\mu \left(\frac{k_x^2}{k_y}-k_y\right)U_0\cos(k_x x)\cos(k_y y){\rm{e}}^{-Qk^2\nu t} \end{split}\right\}$$ (60) 式中, $ U_0 $代表初始速度的幅值; $ k_x=2\text{π}/L_x $ 和 $ k_y=2\text{π}/L_y $分别代表$ x $ 和 $ y $方向上的波数($ k^2=k_x^2+k_y^2 $); $ P_0 $代表参考压力(可以为任意值). 为了简单, 我们采用$k_x = k_y$, $ L_x = L_y=L=1 $以及 $ P_0 =1 $. 此外, $ Q $代表衰减因子, 可以通过调节参数$ Q $来控制流动的衰减速度.
(1) $ Q=1 $时, 代表没有外力, 各项宏观量按标准的泰勒−格林涡流动衰减;
(2) $ Q=0.5 $时, 代表各项宏观量的指数衰减率仅为体力不存在时的一半;
(3) $ Q=0 $时, 体力输入的能量与黏性耗散平衡, 流场中各项宏观量与初始流场保持一致且不随时间变化;
(4) $ Q=-0.5 $时, 体力输入的能量大于黏性耗散, 流场中的压力和速度等宏观量的绝对值将随着时间的增加而增加.
使用$ t=0 $时的解作为流场各宏观量的初始值, 并根据式 (8) 对平衡态分布函数$ g_\alpha^{eq} $进行初始化; 且由于该流动问题对应的初始流场中包含非平衡态部分, 故需要根据式 (27) 和式 (42) 分别对BGK模型和MRT模型下的离散非平衡态分布函数$ g_\alpha^{neq} $初始化; 而两者之和便可以作为初始流场对应的离散分布函数变量$ g_\alpha $.
粗细网格上均采用D2Q9离散速度集. 初始速度幅值$ U_0 $对应的Mach数为$ 0.01 ;$ Reynolds数 $ Re=U_0 L/\nu=100 ;$ 网格数: $ N_x=N_y=100 $. 加密区域范围为$ (x/L,y/L)\in(0.4:0.6,0.1:0.4) $, 加密比例为2. MRT模型下各松弛系数如下
$$ s_e^c=\frac{1}{2\tau_\nu^c},\; s_\epsilon^c=1.7,\; s_q^c=1.9,\; s_\nu^c=\frac{1}{\tau_\nu^c} $$ (62) 在细网格中, MRT模型下的各松弛系数直接使用式 (57) 及 式 (58) 转换而来. 此外还需注意的是原来在计算速度$ {\boldsymbol{u}} $时, 有如下公式(考虑源项为体力项)
$$ \rho u_i =a_{g,i}^{(1)}+\frac{\Delta t}{2}\rho b_i $$ (63) 但本算例中$ {\boldsymbol{b}}=k^2\nu(1-Q){\boldsymbol{u}} $, 其大小与$ {\boldsymbol{u}} $直接相关, 故需要将$ {\boldsymbol{b}} $代入到式 (63) 中来显式得到$ {\boldsymbol{u}} $的计算公式
$$ \rho u_i =a_{g,i}^{(1)}\Bigg/\left[1-\frac{\Delta t}{2} k^2\nu(1-Q)\right] $$ (64) 为了充分验证粗细网格间数据转换关系的正确性, 首先分别绘制了BGK模型和MRT模型在$ Q=0.5 $以及$ tU_0/L=1 $时, 相对压力$ P-P_0 $、水平方向速度$ u_x $以及黏性应力$ \sigma_{xx} $的分布云图. 数值计算结果如图3 所示, 图中的黑色方框为网格加密的界面, 其所包围的区域为加密区域. 从图中可以看出, 各宏观量在加密界面两侧均光滑无间断.
![]()
图 3 强迫泰勒−格林涡各宏观量的分布云图($Q=0.5,\;tU_0/L=1$, 黑色矩形框内的区域为加密区域)Figure 3. The contours of each macroscopic quantity of the forced Taylor-Green vortex ($Q=0.5,\;tU_0/L=1$, the area inside the black rectangular box is the refinement region)此外, 还分别计算了当$ Q $为 $ 0.5 $, $ 0.0 $和$ -0.5 $时, 在$ (x/L,y/L)=(0.4,0.25) $点(粗细网格重合点)处, 相对压力$ P-P_0 $、水平方向速度$ u_x $以及黏性应力$ \sigma_{xx} $随时间的变化, 其数值结果如图4 所示. 在每张子图中, 绘制了在$ Q $取不同的值时, BGK模型以及MRT模型下数值结果与理论解果的对比. 其中实线代表由式 (61) 给出的理论解, 圆圈和倒三角分别代表BGK模型和MRT模型下的数值结果, 而不同的颜色代表$ Q $取不同的值. 观察图4 后发现, 两种碰撞模型的数值结果都与理论解吻合得非常好.
![]()
图 4 强迫泰勒−格林涡中在$ (x/L,y/L)=(0.4,0.25) $点处各宏观量随时间的变化. (a) ~ (c)分别为BGK模型和MRT模型下的相对压力 $ P-P_0 $、水平方向速度 $ u_x $和黏性应力 $ \sigma_{xx} $与理论结果的对比Figure 4. The variation of each macroscopic quantity with time at the point $ (x/L,y/L)=(0.4,0.25) $ in the forced Taylor-Green vortex. (a) ~ (c) are comparisons between the relative pressure $ P-P_0 $, velocity in the $ x $ direction $ u_x $ and viscous stress $ \sigma_{xx} $ under BGK model and MRT model, respectively, and the theoretical results最后, 为了进一步观察粗细网格界面两侧宏观量的连续性, 我们还绘制了当$ Q $为 $ 0.5 $, $ 0.0 $和$ -0.5 $以及$ tU_0/L=1 $时, 在$ y/L=0.25 $处BGK模型和MRT模型下各宏观量随$ x/L $变化的剖线图, 同时也与理论解进行了对比. 其数值结果如图5 所示, 图中两条黑色虚直线间的区域为加密区域. 从图中可以看出, 两种碰撞模型都取得了非常好的数值结果, 且各宏观量在加密界面两侧光滑连续.
BGK模型和MRT模型下良好的数值结果也印证了第 2 节中粗细网格间数据转换关系的正确性. 接下来将以平板泊肃叶流中对流−扩散问题为验证算例, 进一步测试BGK模型下粗细网格间数据转换关系对复杂源项的适应性.
![]()
图 5 强迫泰勒−格林涡中各宏观量的剖线图($tU_0/L=1,\;y/L=0.25$, 两条黑色虚线间的区域为加密区域). (a) ~ (c)分别为BGK模型和MRT模型下相对压力 $ P-P_0 $、水平方向速度 $ u_x $和黏性应力 $ \sigma_{xx} $与理论结果的对比Figure 5. Profile diagram of each macroscopic quantity of forced Taylor-Green vortex ($ tU_0/L=1,\;y/L=0.25 $, the area between the two black dashed lines is refined). (a) ~ (c) are comparisons between the relative pressure $ P-P_0 $, velocity in the $ x $ direction $ u_x $ and viscous stress $ \sigma_{xx} $ under BGK model and MRT model, respectively, and the theoretical results3.2 平板泊肃叶流中对流−扩散问题
如前所述, 合理设计LBM中的源项$ S $可以极大地扩展该方法的适用范围, 许多学者通过CE分析反设计源项$ S $的形式, 以使得修改后的LBM可以满足特定要求[11]. 甚至在很多情况下, LBM被当作一款求解一些特殊偏微分方程的高效通用求解器[42]. 其中, 对流−扩散方程作为一种可以描述由于对流和扩散过程而引起的传递热量、质量或其他物理量输运现象[43]的特殊偏微分方程, 被广泛应用于求解各种对流−扩散问题. 该方程可以被表示如下
$$ \frac{\partial \phi}{\partial t}+\nabla\cdot(\phi {\boldsymbol{u}})=\nabla\cdot(D\nabla\phi) $$ (65) 式中, $ \phi({\boldsymbol{x}},t) $代表一个标量, $ {\boldsymbol{u}}({\boldsymbol{x}},t) $是由N-S方程控制的速度, $ D $为扩散系数.
对于如何使用LBM求解对流−扩散方程, 之前的学者们利用CE分析给出了许多不同的方案[44-48]. 这里使用Chai等[48]提出的通过CE分析反设计源项进而恢复对流−扩散方程的方案, 该方法的标量场中分布函数的演化方程如下
$$ \begin{split} &\phi_{\alpha}\left({\boldsymbol{x}}+{\boldsymbol{\xi}}_{\alpha} \Delta t, t+\Delta t\right)=\left(1-\frac{1}{\tau_\phi}\right)\phi_{\alpha}({\boldsymbol{x}}, t)+\\ &\qquad \frac{1}{\tau_\phi}\phi_\alpha^{eq}({\boldsymbol{x}}, t)+R_\alpha({\boldsymbol{x}}, t) \end{split}$$ (66) 这里$ \phi_\alpha $对应标量场中的分布函数, 而$ R_\alpha $则是为了在恢复式 (65) 时减小额外计算误差所添加的源项, 且$ \phi_\alpha^{eq} $和$ R_\alpha $分别表示如下
$$ \phi_\alpha^{eq}=w_\alpha\phi\left\{1+{\boldsymbol{\xi}}_\alpha\cdot{\boldsymbol{u}}+\frac{1}{2}\left[({\boldsymbol{\xi}}_\alpha\cdot{\boldsymbol{u}})^2-u^2\right]\right\}+\lambda_\alpha\frac{\phi \rho}{\rho_0} $$ (67) $$ R_\alpha=\Delta t\left(1-\frac{1}{2\tau_\phi}\right) w_\alpha{\boldsymbol{\xi}}_\alpha\cdot (\rho\nabla\phi/\rho_0+\phi {\boldsymbol{b}}) $$ (68) 这里$ \nabla\phi $可以直接使用如下公式计算
$$ \nabla\phi=-\frac{2\displaystyle\sum _{\alpha=1}^{q}{\boldsymbol{\xi}}_\alpha(\phi_\alpha-\phi_\alpha^{eq})+\Delta t \phi {\boldsymbol{b}}}{\Delta t (\rho/\rho_0+2\tau_\phi)} $$ (69) 其中$ \tau_\phi=D/\Delta t + 1/2 $以及$\lambda_1=-\displaystyle\sum _{\alpha\neq 1}\lambda_\alpha,\lambda_\alpha=w_\alpha(\alpha\neq 1)$.
接下来将用上述方案去求解平板泊肃叶流中对流−扩散问题. 在计算这个问题时, 需要两套分布函数: $ g_\alpha $用于计算速度场中的密度$ \rho $以及速度$ {\boldsymbol{u}} ; $ $ \phi_\alpha $用于计算标量场中的$ \phi $.
该问题对应的计算简图如图6 所示. 其中上下板处的灰色区域为加密区域, 加密范围为: $ y/H\in\{(0.0:0.1)\cup(0.9:1.0)\} $, 加密比例为$ 2 $. 其他参数为: 定常时中心速度$ U_0 $对应的Mach数为$ 0.01 $; 下板标量: $ \phi_b=0.0 $, 上板标量: $ \phi_t=0.1; $ Reynolds数 $ Re=HU_0/\nu=10, $ $ H $为槽道宽度; Peclet 数 $Pe= HU_0/D=10;$ 网格数 $ N_x=6,N_y=20 $, 粗细网格均采用D2Q9速度集以及BGK碰撞模型. 该流动采用水平方向的体力驱动, 单位质量所受的体力为: $b_x= 8\nu U_0/H^2$. 左右为周期边界条件, 上下边界均使用非平衡态外推边界条件[49], 其中边界点处密度的计算使用非平衡态反弹边界条件[50]中密度的计算方法.
![]()
图 6 平板泊肃叶流中对流−扩散问题示意图Figure 6. Schematic diagram of diffusion problem in a planar Poiseuille flow该流动问题中速度$ u_x(y,t) $所满足的控制方程为
$$ \frac{\partial u_x(y,t)}{\partial t}=b_x+\nu\frac{\partial^2 u_x(y,t)}{\partial^2 y} $$ (70) 对应的解析解可以通过分离变量法求得, 其级数形式可以表示为
$$ \frac{u_x(y,t)}{U_0} =4\left(\frac{y}{H}-\frac{y^2}{H^2}\right)-\sum\limits_{n=0}^{\infty} \frac{32}{(k \text{π})^3} \sin \left(\frac{k \text{π} y}{H}\right) \exp \left(-\frac{k^2 \text{π}^2 \nu t}{H^2}\right) $$ (71) 式中$ k=2 n+1 $. 此外, 标量$ \phi(y,t) $所满足的控制方程可以表示如下
$$ \frac{\partial \phi(y,t)}{\partial t}=D\frac{\partial^2 \phi(y,t)}{\partial y^2} $$ (72) 其对应的解析解也同样可以通过分离变量法给出, 且级数形式表示如下
$$ \begin{split} &\phi(y,t)=\phi_b+(\phi_t-\phi_b)\frac{y}{H}-\\ &\quad \frac{2}{\text{π}}\sum\limits_{k=1}^{\infty}\frac{1}{k}\left[\phi_b+(-1)^{1+k}\phi_t\right]\sin\left(\frac{k \text{π} y}{H}\right)\exp\left(-\frac{k^2\text{π}^2 D t}{H^2}\right)\end{split}\tag{73a} $$ $$ J_x(y,t)=\phi(y,t) u_x(y,t)\tag{73b} $$ $$\begin{split} &J_y(y,t)=-D\Bigg\{\frac{1}{H}(\phi_t-\phi_b)-\\ &\quad \frac{2}{\text{π}}\sum\limits_{k=1}^{\infty}\frac{\text{π}}{H}\left[\phi_b+(-1)^{1+k}\phi_t\right]\cos\left(\frac{k \text{π} y}{H}\right)\exp\left(-\frac{k^2\text{π}^2 D t}{H^2}\right)\Bigg\}\end{split} \tag{73c} $$ 其中$ {\boldsymbol{J}}=-D\nabla\phi+\phi{\boldsymbol{u}} $代表扩散通量与对流通量之和, 而$ J_x $, $ J_y $分别为$ x $, $ y $方向上的分量.
我们使用局部网格加密技术, 分别计算了在$ t^*=t\nu/H^2=0.04,\ 0.1,\ 1.0 $时, 水平方向的速度 $ u_x $、标量 $ \phi $、通量 $ J_x $和$ J_y $并与通过式 (71)和式(73) 得到的理论结果进行了对比. 数值结果如图7 所示, 图中的黑色虚线代表粗细网格的界面, 实线代表理论解, 圆圈代表数值解, 不同的颜色代表不同的时刻. 从图中可以看出, 不同时刻下的数值结果都与理论解吻合得非常好, 而且在加密界面的两侧各宏观量均光滑且连续. 良好的数值结果也进一步验证了第 2 节中粗细网格间数据转换关系的正确性, 同时也验证了该转换关系对复杂源项的良好适应性.
![]()
图 7 平板泊肃叶流中对流−扩散问题(黑色虚线为加密界面). (a) ~ (d)分别为BGK模型下水平方向速度 $ u_x $、标量 $ \phi $、通量 $ J_x $和$ J_y $与理论结果的对比Figure 7. Diffusion problem in a planar Poiseuille flow (the black dashed lines are the refinement interface). (a) ~ (d) represent the velocity in x direction $ u_x $, scalar $ \phi $, flux $ J_y $ and $ J_x $ compared with theoretical results under the BGK model3.3 顶盖驱动方腔流动
顶盖驱动方腔流作为经典的基准算例已被广泛用作数值方法的测试案例. 为了对数值结果进行合理的评价, 使用Tamer等[51]的数值结果作为参考数据并与之进行对比. 使用LBM算法去计算顶盖驱动方腔流时, 如果使用的边界条件较为粗糙且网格分辨率较低时, 顶盖左右两个角点处往往会产生明显的压力“噪声”.
下面使用局部加密算法对该流动进行数值模拟, 用以进一步展示局部加密算法的优势. 本次模拟使用BGK模型以及D2Q9速度集, 雷诺数$Re = U_wL/\nu=1000$, 其中$ L $代表顶盖长度, $ U_w $为顶盖速度其对应的马赫数为$ 0.1 $, $ \nu $为运动黏度系数. 边界条件: 左右壁面以及下壁面均使用修正反弹格式, 顶盖处使用运动边界对应的反弹格式 [52], 即
$$ g_{\bar{\alpha}}({\boldsymbol{x}}_b,t+\Delta t)= g_{\alpha}({\boldsymbol{x}}_b,t+\Delta t)-2w_\alpha\rho_w({\boldsymbol{\xi}}_\alpha\cdot{\boldsymbol{u}}_w) $$ (74) 这里$ g_{\bar{\alpha}} $(流向流场方向的分布函数)代表与$ g_{\alpha} $方向相反的分布函数. 其中$ \rho_w $与$ {\boldsymbol{u}}_w $分别代表边界点处的密度和速度. 其中$ \rho_w $仍采用非平衡态反弹格式[50]中密度的计算方法. 4个角点的边界条件需要额外考虑. 本次模拟中我们对4个角点切线方向的分布函数使用一阶泰勒展开, 偏导数部分使用中心插值, 以右上角点(图8(b))为例
![]()
图 8 顶盖驱动方腔流4个角点示意图, 其中黑色粗实线代表壁面Figure 8. Diagram of four corner points of the square cavity flow driven by the top lid, where the thick black line represents the wall surface$$ \begin{split} &g_{7,9}(x_b,y_b,t)=g_{7,9}(x_b-\Delta x,y_b-\Delta x,t)+\\ &\qquad\frac{1}{2}\left[g_{7,9}(x_b,y_b-\Delta x,t)-g_{7,9}(x_b-2\Delta x,y_b-\Delta x,t)\right]+\\ &\qquad\frac{1}{2}\left[g_{7,9}(x_b-\Delta x,y_b,t)-g_{7,9}(x_b-\Delta x,y_b-2\Delta x,t)\right] \end{split} $$ (75) 其中$ (x_b,y_b) $代表右上角点对应的坐标, 而$ 2,3,6 $方向均使用静止边界的反弹格式. 左上角点也是同样处理.
对整个顶盖附近实施加密, 加密范围为: $ y/L\in(0.7:1.0) $, 加密比例为$ 2 $. 不同分辨率下压力分布云图如图9 所示. 其中, 不仅对比了使用局部加密(UCG-L)和使用均匀粗网格(分辨率与局部加密中的粗网格分辨率一致, UCG)以及均匀细网格(分辨率与局部加密中的细网格分辨率一致, UFG)下的结果, 同时还添加了与UCG-L达到稳态时理论计算耗时相同的均匀网格(EQ)下的结果, 用来对比使用局部网格加密算法的非均匀网格与相同计算量下的均匀网格两者的数值结果. 从图9 中不难看出, 采用局部加密后原来的“噪声”得到了极大的消除, 而且其对“噪声”的抑制效果比相同计算耗时下的均匀网格更好. 这也侧面验证了局部网格加密方法在消除局部“噪声”的有效性. 值得一提的是, Dong等[35]通过详细的分析指出, 这里的压力噪声主要是边界条件隐藏误差导致的, 且可以通过利用双松弛时间碰撞模型优化加以有效消除, 其文章中的数值结果也证明了这一点.
![]()
图 9 $ Re=1000 $时顶盖驱动方腔流的压力分布云图对比, (a) ~ (d)分别为UCG, UFG, EQ和UCG-L对应的压力分布云图 ((d)中的黑色直线为加密界面, 界面以上为加密区域)Figure 9. When $ Re=1000 $, the contours of pressure of the lid-driven cavity flow are compared. (a) ~ (d) are the contours of pressure corresponding to UCG, UFG, EQ and UCG-L, respectively (the black line in (d) is the grid refinement interface, and the area above the black line is the refinement area)图10 中绘制了使用不同分辨率时沿垂直方向穿过几何中心的水平速度$ u_x $随$ y/L $的变化关系以及沿水平方向穿过几何中心的垂直速度$ u_y $随$ x/L $的变化关系. 而且从图10(b) 以及图10(d) 中可以很清楚地看到即使是远离加密区域UCG-L的数值结果仍要好于与其所需计算量相同的均匀网格下的结果, 而且在使用更少计算资源的情况下还取得了与UFG非常相近的数值结果, 这就表明在局部流域执行网格加密可以很好地减小局部数值误差对全流域的影响.
![]()
图 10 $ Re=1000 $时顶盖驱动方腔流中无量纲速度的剖线图, UCG, UFG, EQ和UCG-L与参考数据的对比. (a) 沿垂直直线穿过几何中心的$ u_x/U_w $, (b)为其局部放大, (c) 沿水平直线穿过几何中心的$ u_y/U_w $, (d)为其局部放大Figure 10. Profile diagram of dimensionless velocity in a lid-driven cavity flow when $ Re=1000 $, comparison of UCG, UFG, EQ and UCG-L with reference data. (a) $ u_x/U_w $ in a vertical straight line through the center of the geometry, (b) local amplification of (a), (c) $ u_y/U_w $ in a horizontal straight line through the center of the geometry, (d) local amplification of (c)在表1 中, 我们比较了3种情况下每 1000 个迭代步的 CPU 实际消耗时间. 如果假设在 UCG 上运行到稳态所需的 CPU 时间是$ t $, 若考虑对流时间尺度以及二维情况下则UFG下运行到相同状态下, 理论上所需的 CPU 时间是$ 8 t $, 而上述加密的UCG-L 只需要 $ 3.4 t $, 表1中的测量结果也基本支持上述论断. 值得注意的是, 在本算例中UCG-L所需的计算时间还不及 UFG 的一半, 但其数值结果却与 UFG 非常相近. 有理由相信在某些特殊算例中通过合理布置加密区域, 局部网格加密可以实现更高的加速比.
表 1 不同分辨率下程序演化1000步时CPU消耗时间对比Table 1. Comparison of CPU consumption time when the program evolves 1000 steps at different resolutionsUCG UCG-L UFG CPU time/s 1.656 5.803 6.782 3.4 一维剪切波问题
为了观察局部网格加密技术对LBM算法在数值耗散方面的影响, 我们选择了一维剪切波问题[53]作为研究算例, 其示意图如图11 所示. 该问题对应的控制方程及初始化流场如下
$$ \frac{\partial u_x(y,t)}{\partial t}=\nu\frac{\partial^2 u_x(y,t)}{\partial^2 y} $$ (76) $$ u_x|_{t=0}(y)=u_0\sin (2\text{π} y/H), \quad y\in [0,H]\tag{77a} $$ $$ u_y|_{t=0} = 0 \tag{77b}$$ 式中, $ u_0 $代表初速流场中速度的幅值. 由于物理黏性的存在, 速度$ u_x $会随着时间逐渐衰减, 其对应的理论解为
$$ u_x(y,t)=u_x|_{t=0}(y)\exp[-\nu(2\text{π}/H)^2t] $$ (78) 针对该问题数值模拟的各参数和设置如下: 上下使用周期边界条件; 粗细网格上均采用BGK碰撞模型和D2Q9速度集; 非平衡态分布函数的初始化仍然使用式 (27). 初始速度幅值$ u_0 $对应的Mach数为$ 0.1 $; Reynolds数 $ Re=Hu_0/\nu=10 $, H为槽道宽度; 网格数 $ N_x=4 $, $ N_y=20 $. 分别计算了加密区域范围为: $ (x/L,y/L)\in(0.0:1.0,0.1:0.25), $ $(x/L,y/L)\in (0.0: 1.0, 0.1:0.4),$ $ (x/L,y/L)\in(0.0:1.0,0.1:0.6) $ 和 $(x/L, y/L)\in (0.0:1.0,0.4:0.6)$, 且加密比例均为$ 2 $的情况下, UCG-L和UCG下的数值耗散$ \nu_N $与物理黏性$ \nu $之比, 其中数值黏性可以通过反解式 (78) 求得
$$ \nu_N(y,t)=\ln\left(\frac{u_x|_{t=0}(y)}{u_x^N(y,t)}\right)\Bigg/\left[(2\text{π}/H)^2t\right] $$ (79) 这里$ u_x^N $代表数值结果.
分别绘制出UCG-L和UCG下的$ \nu_N/\nu $随着$ y/H $变化的剖线图, 其数值结果如图12 所示. 图12 展示了在不同区域使用局部网格加密后, 在给定时刻($t^*=t\nu/H^2=4\text{π}^2/\ln 2$)下的数值黏性与物理黏性之比(红色线)以及均匀粗网格下的数值黏性与物理黏性之比(蓝色线), 图中的黑色虚线代表粗细网格的分界线. 从图中可以发现, 不同加密区域下, 加密界面两侧的数值黏性相比于均匀网格下, 表现出了不一样的结果, 既有增大的情况也有降低的情况. 而其中数值黏性小于物理黏性现象的出现可能与加密界面附近的非物理震荡[19,54-55]有着密切的关系, 故还需要一系列更加系统的理论研究将局部网格加密技术引起的数值耗散和数值色散与加密区域附近的非物理震荡联系起来, 而这也将作为我们后续工作的重点.
![]()
图 12 一维剪切波问题在不同加密区域下均匀网格和非均匀网格对应的数值黏性和物理黏性之比的分布情况(黑色虚线为加密界面). (a) ~ (d)分别为不同加密区域下UCG-L和UCG对应的数值黏性和物理黏性之比随$ y/H $变化的剖线图 (续)Figure 12. For one-dimensional shear wave problem, the distribution of the ratio of numerical viscosity and physical viscosity corresponding to uniform and non-uniform grids in different refinement regions (the black dashed line is the refinement interface). (a) ~ (d) are the cross-section plots of the ratio of numerical viscosity and physical viscosity corresponding to UCG-L and UCG with $ y/H $ in different refinement regions, respectively (continued)4. 总结
本文在考虑任意源项存在的前提下, 提出了一套规范且简洁的粗细网格在重合点处的数据转换关系推导方法, 该方法既可以用于BGK模型也可以适用于MRT模型, 同时也容易被推广到其他碰撞模型下的数据转换. 该推导方法从保证粗细网格重合点处对应的连续分布函数一致出发, 除忽略离散误差外, 并不依赖于CE展开以及非平衡态近似. 这也说明, 局部网格加密算法不仅适用于连续流, 而且也适用于高阶非平衡态流动问题. 此外, 我们还从理论上证明了保证粗细网格间非平衡态部分的一阶CE近似一致, 便可以保证整个非平衡态部分的一致, 这也给之前学者对非平衡态部分做一阶CE近似提供了理论依据. 同时, 我们还通过引入新的Hermite转换矩阵$ {\boldsymbol{M}}_H $直接构建了MRT模型下, 时空离散前后分布函数对应的各阶Hermite矩间的关系, 从而将BGK模型下基于Hermite展开的初始化方法推广到了MRT模型下.
为了验证包含复杂源项的粗细网格间数据转换关系, 考虑了两个具体算例. 第1个算例是外力作用下的不可压强迫泰勒−格林涡流动, 这个问题涉及非均匀非定常流场, 并且由源项表示的外力也为非均匀场. 初始流场对应的分布函数采用Hermite展开来准确构建. 流场内部的部分区域采用细网格. 计算结果显示粗细网格边界处压力、速度、应力场均光滑衔接, 这些宏观量随时间和空间的演化剖面都与解析解完全吻合, 并且BGK碰撞模型和MRT碰撞模型下的计算结果也基本完全一致, 证明了我们推导的基于两种碰撞模型转换关系的正确性. 第2个算例是平板泊肃叶流驱动的非定常对流−扩散问题, 这里标量场对应的Boltzmann方程需要引进源项, 才能准确恢复宏观非定常对流−扩散方程. 在上下两个壁面附近, 我们采用细网格局部加密. 本问题的标量场在流向有对流通量, 在垂直壁面方向有扩散通量, 并且它们分布都不均匀, 且随时间演化. 通过分离变量法, 给出了标量场的非定常解析解. 计算结果显示速度场、标量场和标量通量的两个分量在粗细网格界面都光滑过渡, 并与解析解完全吻合, 进一步验证了一般源项下转换关系的正确性. 此外, 我们还使用局部网格加密技术对经典的顶盖驱动方腔流动问题进行了数值模拟, 数值结果表明局部网格加密技术在处理局部压力噪声方面有着明显的优势, 而且通过对比使用局部加密的非均匀网格和不使用局部加密的均匀网格的计算耗时以及宏观速度的剖线图, 可以较为清楚地看到局部网格加密技术在处理复杂流动问题时的优势. 最后, 为了观察局部网格加密技术对LBM算法在数值耗散方面的影响, 我们通过在一维剪切波中的不同区域使用局部网格加密, 并将其数值黏性与均匀粗网格下的数值黏性进行了对比, 发现由局部加密引起的数值黏性与加密区域的选取有很大的关系, 且在部分加密区域下还出现了数值黏性小于物理黏性的情况, 一系列更加精细的算例和理论分析需要被考虑来阐述该现象的发生.
在上述的数值验证中, 我们的算例均为二维算例, 并没有涉及更加复杂的三维流动问题, 这主要是由于LBM中的局部网格加密技术经过20多年的发展已经趋于成熟, 相关的文献中 [25,37,56]均已报道了大量使用局部网格加密技术计算复杂三维流动问题的测试结果, 相关的数值结果均展示了局部网格加密技术在三维复杂流动问题中的优势及适应性. 本文的主要创新点在于从粗细网格间分布函数需要转换的本质出发, 构建了一套直观且简洁的包含任意源项下分布函数转换关系的推导过程, 且我们的推导结果与前人的结果完全一致, 这一点在 2.1 中也给出了解释. 而如上所述, LBM中的局部网格加密技术已经被广泛应用到各种复杂流动问题的计算中, 且其使用的分布函数转换关系与我们的也完全一致, 故本文中并没有对更加复杂的三维流动进行模拟, 而是借用二维情况下的强迫泰勒−格林涡流动和平板泊肃叶流中的对流−扩散问题, 着重考察了分布函数转换关系对复杂源项的适应性, 以及通过对顶盖驱动方腔流动的模拟初步展示局部网格加密技术在处理复杂流动问题方面的优势.
同时需要注意的是, 粗细网格重合点处的分布函数转换关系式 (48)和式(49) 中的粗细网格分辨率之比 $ n $ 理论上可以采用任意正整数, 但本文粗细网格间分布函数的转换关系是在忽略时空离散误差的前提下得出的, 而实际情况下由于离散误差的存在, 不同网格分辨率下得到的数值结果中所包含的离散误差也是不一样的, 故粗细网格界面处的数值结果间往往会产生一定的“数值间断”, 且该“数值间断”随着加密比例的增加而增加, 对于一些较为复杂的流动, 这往往会在粗细网格界面产生非物理的震荡, 故上面的算例中局部网格的加密倍数均为$ 2 $, 而且绝大多数文献的网格加密比例在计算中往往也被限定为$ 2 $.
必须指出的是, 本文的分析并没有考虑粗细网格间不同数值离散误差有可能在网格加密界面处产生的数值震荡. 这个问题在已有文献中有一些讨论[19,57], 但需要进一步做理论分析并找到消除数值震荡的有效方法. 其中, 不同碰撞模型对粗细网格界面的数值稳定性有着不同的影响[19], 但也需要从理论上进一步做分析. 如果网格加密边界与物理壁面相交, 那么物理壁面的反弹格式局部误差[35]与粗细网格界面误差交汇, 此时往往会在壁面附近产生一定的数值震荡[54], 该问题的理论分析会更具挑战性, 但对实际应用问题却非常重要. 这些都是需要进一步研究的问题.
此外, 测试和分析局部网格加密算法在数值色散、数值耗散方面的表现对于该技术的实际应用有着非常重要的意义, 尤其是有助于加密界面处非物理震荡问题[19,54-55]的解决. 在LBM的局部网格加密算法中, 粗细网格上仍然使用标准的LBM算法, 而不同的是需要在粗细网格交界面处进行合理的数据转换. 其中, 数据转换时往往涉及未知分布函数的构建, 目前构建方式可以分为有限差分方法和有限体积方法两类, 而不同的构建方式对局部网格加密算法的数值耗散和数值色散都有不同的影响. 故研究局部网格加密算法对LBM在数值色散和数值耗散方面的影响, 不仅需要详细对比不同的构建格式在数值耗散和数值色散方面的表现, 还需要从理论上分析不同构建格式间结果差异的原因, 以及数值耗散和数值色散对加密界面处非物理震荡的影响, 而这也将是我们下一步的工作重点. 考虑到该问题的复杂性, 故我们将在后续文章中尝试对该问题做较为完整和详细的介绍.
附录A MRT模型下时空离散前后的分布函数对应的各阶Hermite间的关系
首先, 需要对式 (29) 两边同时左乘Hermite多项式矩阵$ {\boldsymbol{H}} $
$$ {\boldsymbol{H}}|f\rangle={\boldsymbol{H}}\left({\boldsymbol{I}}+\frac{\Delta t}{2}{\boldsymbol{\varLambda}}\right)^{-1}{\boldsymbol{H}}^{-1}\left({\boldsymbol{H}}|g^{neq}\rangle+\frac{\Delta t}{2} {\boldsymbol{H}}|S\rangle\right)+{\boldsymbol{H}}|g^{eq}\rangle \tag{A1} $$ 为了让上式可以直接反映连续分布函数变量$ |f\rangle $的各阶Hermite矩与离散分布函数变量$ |g\rangle $各阶Hermite矩间的关系, 上式中的Hermite多项式矩阵$ {\boldsymbol{H}} $可以直接由Hermite多项式组成, 以D2Q9速度集为例
$$ \begin{array}{*{20}{l}} {\boldsymbol{H}}=\left[\begin{array}{c} \langle H^{(0)}|\\ \langle H^{(1)}_x|\\ \langle H^{(1)}_y|\\ \langle H^{(2)}_{xx}|\\ \langle H^{(2)}_{xy}|\\ \langle H^{(2)}_{yy}|\\ \langle H^{(3)}_{xxy}|\\ \langle H^{(3)}_{xyy}|\\ \langle H^{(4)}_{xxyy}|\\ \end{array}\right]=\left[\begin{array}{ccccccccc} 1&1&1&1&1&1&1&1&1 \\ 0&\sqrt{3}&0&-\sqrt{3}&0&\sqrt{3}&-\sqrt{3}&-\sqrt{3}&\sqrt{3} \\ 0&0&\sqrt{3}&0&-\sqrt{3}&\sqrt{3}&\sqrt{3}&-\sqrt{3}&-\sqrt{3} \\ -1&2&-1&2&-1&2&2&2&2 \\ 0&0&0&0&0&3&-3&3&-3 \\ -1&-1&2&-1&2&2&2&2&2 \\ 0&0&-\sqrt{3}&0&\sqrt{3}&2\sqrt{3}&2\sqrt{3}&-2\sqrt{3}&-2\sqrt{3} \\ 0&-\sqrt{3}&0&\sqrt{3}&0&2\sqrt{3}&-2\sqrt{3}&-2\sqrt{3}&2\sqrt{3} \\ 1&-2&-2&-2&-2&4&4&4&4 \end{array}\right] \end{array} \tag{A2}$$ 上式中各阶Hermite 多项式的具体形式可以表示如下[58]
$$ \left. \begin{split} &H^{(0)}=1\\ & H^{(1)}_x({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,x}\\ &H^{(1)}_y({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,y}\\ &H^{(2)}_{xx}({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,x}\xi_{\alpha,x}-1\\ &H^{(2)}_{xy}({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,x}\xi_{\alpha,y} \\ &H^{(2)}_{yy}({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,y}\xi_{\alpha,y}-1\\ &H^{(3)}_{xxy}({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,x}^2\xi_{\alpha,y}-\xi_{\alpha,y} \\ &H^{(3)}_{xyy}({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,x}\xi_{\alpha,y}^2-\xi_{\alpha,x}\\ &H^{(4)}_{xxyy}({\boldsymbol{\xi}}_\alpha)=\xi_{\alpha,x}^2\xi_{\alpha,y}^2-(\xi_{\alpha,x}^2+\xi_{\alpha,y}^2)+1\end{split} \right\} \tag{A3} $$ 此时, 式 (A1) 中的$ {\boldsymbol{H}}|f\rangle $可以表示如下
$$ \begin{array}{*{20}{l}} {\boldsymbol{H}}|f\rangle=\left[\begin{array}{c} \displaystyle\sum _{\alpha=1}^{q} f_\alpha H^{(0)({\boldsymbol{\xi}}_\alpha)}\\ \displaystyle\sum _{\alpha=1}^{q} f_\alpha H^{(1)}_x({\boldsymbol{\xi}}_\alpha)\\ \displaystyle\sum _{\alpha=1}^{q} f_\alpha H^{(1)}_y({\boldsymbol{\xi}}_\alpha)\\ \displaystyle\sum _{\alpha=1}^{q} f_\alpha H^{(2)}_{xx}({\boldsymbol{\xi}}_\alpha)\\ \displaystyle\sum _{\alpha=1}^{q} f_\alpha H^{(2)}_{xy}({\boldsymbol{\xi}}_\alpha)\\ \displaystyle\sum _{\alpha=1}^{q} f_\alpha H^{(2)}_{yy}({\boldsymbol{\xi}}_\alpha)\\ \displaystyle\sum _{\alpha=1}^{q} f_\alpha H^{(3)}_{xxy}({\boldsymbol{\xi}}_\alpha)\\ \displaystyle\sum _{\alpha=1}^{q} f_\alpha H^{(3)}_{xyy}({\boldsymbol{\xi}}_\alpha)\\ \displaystyle\sum _{\alpha=1}^{q} f_\alpha H^{(4)}_{xxyy}({\boldsymbol{\xi}}_\alpha)\\ \end{array}\right]=\left[\begin{array}{c} a^{(0)}\\ a^{(1)}_x\\ a^{(1)}_y\\ a^{(2)}_{xx}\\ a^{(2)}_{xy}\\ a^{(2)}_{yy}\\ a^{(3)}_{xxy}\\ a^{(3)}_{xyy}\\ a^{(4)}_{xxyy}\\ \end{array}\right]=|{{a}}\rangle \end{array} \tag{A4}$$ 同理, 式 (A1) 中的${\boldsymbol{H}}|g^{neq}\rangle = |{{a}}_g^{neq}\rangle$, ${\boldsymbol{H}}|S\rangle = |{{a}}_s\rangle$和${\boldsymbol{H}}|g^{eq}\rangle = |{{a}}_g^{eq}\rangle$可以分别表示如下
$$ \begin{array}{*{20}{l}} |{{a}}_g^{neq}\rangle=\left[\begin{array}{c} {{a}}^{neq,(0)}_{g}\\ a^{neq,(1)}_{g,x}\\ a^{neq,(1)}_{g,y}\\ a^{neq,(2)}_{g,xx}\\ a^{neq,(2)}_{g,xy}\\ a^{neq,(2)}_{g,yy}\\ a^{neq,(3)}_{g,xxy}\\ a^{neq,(3)}_{g,xyy}\\ a^{neq,(4)}_{g,xxyy}\\ \end{array}\right] ,\;\; |{{a}}_s\rangle=\left[\begin{array}{c} a^{(0)}_{s}\\ a^{(1)}_{s,x}\\ a^{(1)}_{s,y}\\ a^{(2)}_{s,xx}\\ a^{(2)}_{s,xy}\\ a^{(2)}_{s,yy}\\ a^{(3)}_{s,xxy}\\ a^{(3)}_{s,xyy}\\ a^{(4)}_{s,xxyy}\\ \end{array}\right] , \;\; |{{a}}_g^{eq}\rangle=\left[\begin{array}{c} a^{eq,(0)}_{g}\\ a^{eq,(1)}_{g,x}\\ a^{eq,(1)}_{g,y}\\ a^{eq,(2)}_{g,xx}\\ a^{eq,(2)}_{g,xy}\\ a^{eq,(2)}_{g,yy}\\ a^{eq,(3)}_{g,xxy}\\ a^{eq,(3)}_{g,xyy}\\ a^{eq,(4)}_{g,xxyy}\\ \end{array}\right] \end{array}\tag{A5} $$ 进而式 (A1) 可以表示为
$$ |{{a}}\rangle={\boldsymbol{H}}\left({\boldsymbol{I}}+\frac{\Delta t}{2}{\boldsymbol{\varLambda}}\right)^{-1}{\boldsymbol{H}}^{-1}\left[|{{a}}_g^{neq}\rangle+\frac{\Delta t}{2} |{{a}}_s\rangle\right]+|{{a}}_g^{eq}\rangle \tag{A6} $$ 上式直接给出了连续分布函数变量$ |f\rangle $和离散分布函数变量$ |g\rangle $对应的Hermite矩间的关系. 在对式 (32) 两边取逆后, 可以得到
$$ {\boldsymbol{\varLambda}}^{-1}=[\Delta t({\boldsymbol{M}}^{-1} {\boldsymbol{\hat{S}}}{\boldsymbol{M}})^{-1}-\frac{\Delta t}{2}{\boldsymbol{I}}] \tag{A7}$$ 若记Hermite转换矩阵$ {\boldsymbol{M}}_H $为如下形式
$$ \left[{\boldsymbol{I}}-\frac{1}{2}({\boldsymbol{M}}{\boldsymbol{H}}^{-1})^{-1}\hat{{\boldsymbol{S}}}({\boldsymbol{M}}{\boldsymbol{H}}^{-1})\right]^{-1}\equiv {\boldsymbol{M}}_H \tag{A8} $$ 结合式 (A7) 与式 (32) 后, 式 (A6) 可以被进一步简化如下
$$ |a\rangle={\boldsymbol{M}}_H^{-1}\left(|a_g^{neq}\rangle+\frac{\Delta t}{2} |a_s\rangle\right)+|a_g^{eq}\rangle \tag{A9}$$ 上式两边同时减去平衡态分布函数对应的各阶Hermite矩后, 便可以得到连续非平衡态分布函数变量$ |f^{neq}\rangle $的各阶Hermite矩$ |a^{neq}\rangle $与离散非平衡态分布函数变量$ |g^{neq}\rangle $的各阶Hermite矩$ |a_g^{neq}\rangle $间的关系如下
$$ |a^{neq}\rangle={\boldsymbol{M}}_H^{-1}\left(|a_g^{neq}\rangle+\frac{\Delta t}{2} |a_s\rangle\right) \tag{A10}$$ D2Q9速度集下的Hermite转换矩阵$ {\boldsymbol{M}}_H $可以显式表示为
$$ {\boldsymbol{M}}_H=\left[\begin{array}{ccccccccc} 1&0&0&0&0&0&0&0&0 \\ 0&1&0&0&0&0&0&0&0 \\ 0&0&1&0&0&0&0&0&0 \\ c_{21}/C_2&0&0&c_{22}/C_2&0&c_{23}/C_2&0&0&0 \\ 0&0&0&0& d_2 &0&0&0&0 \\ c_{21}/C_2&0&0&c_{23}/C_2&0&c_{22}/C_2&0&0&0 \\ 0&0&c_{31}&0&0&0&d_{3}&0&0 \\ 0&c_{31}&0&0&0&0&0&d_{3}&0 \\ c_{41}/C_2&0&0&c_{42}/C_2&0&c_{42}/C_2&0&0&d_4 \end{array}\right] \tag{A11}$$ 上述Hermite转换矩阵$ {\boldsymbol{M}}_H $中各参数表示的具体表达式如下
$$ \left. \begin{aligned} &C_2=\left(1-\frac{s_\epsilon}{2}\right)\left(1-\frac{s_e}{2}-\frac{s_\nu}{2}+\frac{s_es_\nu}{4}\right)\\ &c_{21}=-\frac{s_e}{2}\left(1-\frac{s_\epsilon}{2}-\frac{s_\nu}{2}+\frac{s_\epsilon s_\nu}{4}\right)\\ &c_{22}=1-\frac{s_e}{4}-\frac{s_\epsilon}{2}+\frac{s_es_\epsilon}{8}-\frac{s_\nu}{4}+\frac{s_\epsilon s_\nu}{8} \\ &c_{23}=\frac{1}{4}\left(s_e-\frac{s_e s_\epsilon}{2}-s_\nu+\frac{s_\epsilon s_\nu}{2}\right)\\ &\ d_2=\frac{2}{2-s_\nu} \end{aligned} \right\} \tag{A12}$$ $$ \left. \begin{aligned} &c_{31}=-\frac{s_q}{2-s_q}\\ &d_{3}=\frac{2}{2-s_q} \end{aligned} \right\} \tag{A13} $$ $$ \left. \begin{aligned} &c_{41}=-s_e+\frac{s_\epsilon}{2}+\frac{s_e s_\epsilon}{4}+\frac{s_es_\nu}{2}-\frac{s_\epsilon s_\nu}{4}-\frac{s_e s_\epsilon s_\nu}{8}\\ &c_{42}=\frac{1}{2}\left(s_e-s_\epsilon-\frac{s_es_\nu}{2}+\frac{s_\epsilon s_\nu}{2}\right)\\ &d_{4}=\frac{2}{2-s_\epsilon} \end{aligned} \right\}\tag{A14} $$ 至此, 可以根据式 (A10) 及式 (A11) 将$ {{a}}_g^{neq} $的各阶Hermite矩显式表示为
$$ \left. \begin{split} &a^{neq,(0)}_g=-\frac{\Delta t}{2}a_s^{(0)}\\ &a^{neq,(1)}_{g,x}=-\frac{\Delta t}{2}a_{s,x}^{(1)}, \ a^{neq,(1)}_{g,y}=-\frac{\Delta t}{2}a_{s,y}^{(1)}\\ &a^{neq,(2)}_{g,xx}=-\frac{1}{C_2}\left(c_{22}\sigma_{xx}+c_{23}\sigma_{yy}\right)-\frac{\Delta t}{2}a_{s,xx}^{(2)}\\ &a^{neq,(2)}_{g,xy}=-d_2 \sigma_{xy}-\frac{\Delta t}{2}a_{s,xy}^{(2)}\\ &a^{neq,(2)}_{g,yy}=-\frac{1}{C_2}\left(c_{23}\sigma_{xx}+c_{22}\sigma_{yy}\right)-\frac{\Delta t}{2}a_{s,yy}^{(2)}\\ &a^{neq,(3)}_{g,xxy}=d_3a^{neq,(3)}_{xxy}-\frac{\Delta t}{2}a_{s,xxy}^{(3)}\\ &a^{neq,(3)}_{g,xyy}=d_3a^{neq,(3)}_{xyy}-\frac{\Delta t}{2}a_{s,xyy}^{(3)}\\ &a^{neq,(4)}_{g,xxyy}=-\frac{c_{42}}{C_2}\left(\sigma_{xx}+\sigma_{yy}\right)+d_4 a^{neq,(4)}_{xxyy}-\frac{\Delta t}{2}a_{s,xxyy}^{(4)} \end{split} \right\} \tag{A15}$$ 当考虑源项为体力项时, $ |g^{neq}\rangle $的前6个矩的具体形式如下
$$ \left. \begin{split} &a^{neq,(0)}_g=0\\ &a^{neq,(1)}_{g,x}=-\frac{\Delta t}{2}\rho b_x\\ &a^{neq,(1)}_{g,y}=-\frac{\Delta t}{2}\rho b_y\\ &a^{neq,(2)}_{g,xx}=-\frac{1}{C_2}\left(c_{22}\sigma_{xx}+c_{23}\sigma_{yy}\right)-\Delta t\rho b_x u_x\\ &a^{neq,(2)}_{g,xy}=-d_2 \sigma_{xy}-\frac{\Delta t}{2}\rho \left(b_x u_y+b_y u_x\right)\\ &a^{neq,(2)}_{g,yy}=-\frac{1}{C_2}\left(c_{23}\sigma_{xx}+c_{22}\sigma_{yy}\right)-\Delta t\rho b_y u_y \end{split} \right\} \tag{A16} $$ 此外, 对比式 (24) 和式 (A16) 后可以发现: 当MRT模型中与体积黏度系数有关的格子松弛参数$ s_e $被设为 $ s_e=s_\nu $时, MRT模型中的密度$ \rho $、速度$ {\boldsymbol{u}} $、黏性应力$ {\boldsymbol{\sigma}} $以及初始化离散分布函数$ g_\alpha $(Hermite展开到两阶)的计算表达式与BGK模型下完全一致.
附录B 各阶连续非平衡态分布函数的递推关系
直接从BGK模型下的DVBE出发
$$ \frac{{\rm{D}}}{{\rm{D}}t}f_\alpha\equiv\partial _t f_\alpha + {\boldsymbol{\xi}}_\alpha\cdot \nabla f_\alpha=\frac{1}{\tau}(f_\alpha^{eq}-f_\alpha)+S_\alpha \tag{B1}$$ 进行简单移项后
$$ f_\alpha=f_\alpha^{eq}-\tau\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha + \tau S_\alpha \tag{B2}$$ 对式 (B2) 使用若干次麦克斯韦迭代
$$ \begin{split} &f_\alpha=f_\alpha^{eq}-\tau\left(\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha^{eq}- S_\alpha\right)+\tau^2\frac{{\rm{D}}}{{\rm{D}}t}\left(\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha-S_\alpha\right)= \\ &\qquad f_\alpha^{eq}-\tau\left(\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha^{eq}- S_\alpha\right)+\tau^2\frac{{\rm{D}}}{{\rm{D}}t}\left(\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha^{eq}-S_\alpha\right)-\\ &\qquad \tau^3\left(\frac{{\rm{D}}}{{\rm{D}}t}\right)^2\left(\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha^{eq}-S_\alpha\right)+O(\tau^4)\end{split} \tag{B3}$$ 上式也可以写作如下形式
$$ f_\alpha^{neq}=f_\alpha-f_\alpha^{eq}=\sum\limits_{k=1}^{N}\tau^k f_\alpha^{(k)}+O(\tau^{N+1}) \tag{B4}$$ 其中$ f_\alpha^{(k)} $的前几阶可以表示如下
$$ f_\alpha^{(1)}=-\left(\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha ^{eq}-S_\alpha\right) \tag{B5a}$$ $$ f_\alpha^{(2)}=\frac{{\rm{D}}}{{\rm{D}}t}\left(\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha ^{eq}-S_\alpha\right)=-\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha^{(1)}\tag{B5b} $$ $$ f_\alpha^{(3)}=-\left(\frac{{\rm{D}}}{{\rm{D}}t}\right)^2\left(\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha ^{eq}-S_\alpha\right)=\left(\frac{{\rm{D}}}{{\rm{D}}t}\right)^2f_\alpha^{(1)}=-\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha^{(2)} \tag{B5c}$$ 归纳可得如下关系
$$ f_\alpha^{(k)}=-\frac{{\rm{D}}}{{\rm{D}}t}f_\alpha^{(k-1)}=(-1)^{k-1}\left(\frac{{\rm{D}}}{{\rm{D}}t}\right)^{k-1}f_\alpha^{(1)} \tag{B6}$$ 从上式不难看出, 各阶$ f_\alpha^{(k)} $间存在明显的递推关系, 且均可以通过$ f_\alpha^{(1)} $进行表示.
-
![]()
图 2 粗细网格排布示意图. 灰色区域为粗细网格的重叠区域
Figure 2. Schematic diagram of the coarse and fine grid arrangement. The grey area is the overlap between the coarse and fine grids
![]()
图 3 强迫泰勒−格林涡各宏观量的分布云图($Q=0.5,\;tU_0/L=1$, 黑色矩形框内的区域为加密区域)
Figure 3. The contours of each macroscopic quantity of the forced Taylor-Green vortex ($Q=0.5,\;tU_0/L=1$, the area inside the black rectangular box is the refinement region)
![]()
图 4 强迫泰勒−格林涡中在$ (x/L,y/L)=(0.4,0.25) $点处各宏观量随时间的变化. (a) ~ (c)分别为BGK模型和MRT模型下的相对压力 $ P-P_0 $、水平方向速度 $ u_x $和黏性应力 $ \sigma_{xx} $与理论结果的对比
Figure 4. The variation of each macroscopic quantity with time at the point $ (x/L,y/L)=(0.4,0.25) $ in the forced Taylor-Green vortex. (a) ~ (c) are comparisons between the relative pressure $ P-P_0 $, velocity in the $ x $ direction $ u_x $ and viscous stress $ \sigma_{xx} $ under BGK model and MRT model, respectively, and the theoretical results
![]()
图 5 强迫泰勒−格林涡中各宏观量的剖线图($tU_0/L=1,\;y/L=0.25$, 两条黑色虚线间的区域为加密区域). (a) ~ (c)分别为BGK模型和MRT模型下相对压力 $ P-P_0 $、水平方向速度 $ u_x $和黏性应力 $ \sigma_{xx} $与理论结果的对比
Figure 5. Profile diagram of each macroscopic quantity of forced Taylor-Green vortex ($ tU_0/L=1,\;y/L=0.25 $, the area between the two black dashed lines is refined). (a) ~ (c) are comparisons between the relative pressure $ P-P_0 $, velocity in the $ x $ direction $ u_x $ and viscous stress $ \sigma_{xx} $ under BGK model and MRT model, respectively, and the theoretical results
![]()
图 6 平板泊肃叶流中对流−扩散问题示意图
Figure 6. Schematic diagram of diffusion problem in a planar Poiseuille flow
![]()
图 7 平板泊肃叶流中对流−扩散问题(黑色虚线为加密界面). (a) ~ (d)分别为BGK模型下水平方向速度 $ u_x $、标量 $ \phi $、通量 $ J_x $和$ J_y $与理论结果的对比
Figure 7. Diffusion problem in a planar Poiseuille flow (the black dashed lines are the refinement interface). (a) ~ (d) represent the velocity in x direction $ u_x $, scalar $ \phi $, flux $ J_y $ and $ J_x $ compared with theoretical results under the BGK model
![]()
图 8 顶盖驱动方腔流4个角点示意图, 其中黑色粗实线代表壁面
Figure 8. Diagram of four corner points of the square cavity flow driven by the top lid, where the thick black line represents the wall surface
![]()
图 9 $ Re=1000 $时顶盖驱动方腔流的压力分布云图对比, (a) ~ (d)分别为UCG, UFG, EQ和UCG-L对应的压力分布云图 ((d)中的黑色直线为加密界面, 界面以上为加密区域)
Figure 9. When $ Re=1000 $, the contours of pressure of the lid-driven cavity flow are compared. (a) ~ (d) are the contours of pressure corresponding to UCG, UFG, EQ and UCG-L, respectively (the black line in (d) is the grid refinement interface, and the area above the black line is the refinement area)
![]()
图 10 $ Re=1000 $时顶盖驱动方腔流中无量纲速度的剖线图, UCG, UFG, EQ和UCG-L与参考数据的对比. (a) 沿垂直直线穿过几何中心的$ u_x/U_w $, (b)为其局部放大, (c) 沿水平直线穿过几何中心的$ u_y/U_w $, (d)为其局部放大
Figure 10. Profile diagram of dimensionless velocity in a lid-driven cavity flow when $ Re=1000 $, comparison of UCG, UFG, EQ and UCG-L with reference data. (a) $ u_x/U_w $ in a vertical straight line through the center of the geometry, (b) local amplification of (a), (c) $ u_y/U_w $ in a horizontal straight line through the center of the geometry, (d) local amplification of (c)
![]()
图 12 一维剪切波问题在不同加密区域下均匀网格和非均匀网格对应的数值黏性和物理黏性之比的分布情况(黑色虚线为加密界面). (a) ~ (d)分别为不同加密区域下UCG-L和UCG对应的数值黏性和物理黏性之比随$ y/H $变化的剖线图 (续)
Figure 12. For one-dimensional shear wave problem, the distribution of the ratio of numerical viscosity and physical viscosity corresponding to uniform and non-uniform grids in different refinement regions (the black dashed line is the refinement interface). (a) ~ (d) are the cross-section plots of the ratio of numerical viscosity and physical viscosity corresponding to UCG-L and UCG with $ y/H $ in different refinement regions, respectively (continued)
表 1 不同分辨率下程序演化1000步时CPU消耗时间对比
Table 1 Comparison of CPU consumption time when the program evolves 1000 steps at different resolutions
UCG UCG-L UFG CPU time/s 1.656 5.803 6.782 
下载: 导出CSV
-
[1] Kerimo J, Girimaji SS. Boltzmann-BGK approach to simulating weakly compressible 3D turbulence: comparison between lattice Boltzmann and gas kinetic methods. Journal of Turbulence, 2007, 8: N46
[2] Lin CD, Xu AG, Zhang GC, et al. Double-distribution-function discrete Boltzmann model for combustion. Combustion and Flame, 2016, 164: 137-151 doi: 10.1016/j.combustflame.2015.11.010
[3] Shan XW, Chen HD. Lattice Boltzmann model for simulating flows with multiple phases and components. Physical Review E, 1993, 47(3): 1815 doi: 10.1103/PhysRevE.47.1815
[4] Dellar PJ. Electromagnetic waves in lattice Boltzmann magnetohydrodynamics. Europhysics Letters, 2010, 90(5): 50002 doi: 10.1209/0295-5075/90/50002
[5] Aidun CK, Clausen JR. Lattice-Boltzmann method for complex flows. Annual Review of Fluid Mechanics, 2010, 42: 439-472 doi: 10.1146/annurev-fluid-121108-145519
[6] Lallemand P, Luo LS, Krafczyk M, et al. The lattice Boltzmann method for nearly incompressible flows. Journal of Computational Physics, 2021, 431: 109713 doi: 10.1016/j.jcp.2020.109713
[7] Marié S, Ricot D, Sagaut P. Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics. Journal of Computational Physics, 2009, 228(4): 1056-1070 doi: 10.1016/j.jcp.2008.10.021
[8] Hui X, Pierre S. Optimal low-dispersion low-dissipation lbm schemes for computational aeroacoustics. Journal of Computational Physics, 2011, 230(13): 5353-5382 doi: 10.1016/j.jcp.2011.03.040
[9] Kuznik F, Obrecht C, Rusaouen G, et al. LBM based flow simulation using gpu computing processor. Computers & Mathematics with Applications, 2010, 59(7): 2380-2392
[10] Blair S, Albing C, Grund A, et al. Accelerating an mpi lattice Boltzmann code using openacc//Proceedings of the Second Workshop on Accelerator Programming Using Directives, 2015: 1-9
[11] Chen T, Wang LP, Lai J, et al. Inverse design of mesoscopic models for compressible flow using the Chapman-Enskog analysis. Advances in Aerodynamics, 2021, 3(1): 1-25 doi: 10.1186/s42774-020-00055-6
[12] 陈涛. 可压缩湍流介观模型设计及应用. [博士论文]. 北京: 北京大学, 2021 Chen Tao. Design and application of mesoscopic models for compressible turbulent flow. [PhD Thesis]. Beijing: Peking University (in Chinese)
[13] Filippova O, Hänel D. Grid refinement for lattice-BGK models. Journal of Computational Physics, 1998, 147(1): 219-228 doi: 10.1006/jcph.1998.6089
[14] Lin CL, Lai YG. Lattice Boltzmann method on composite grids. Physical Review E, 2000, 62(2): 2219 doi: 10.1103/PhysRevE.62.2219
[15] Dupuis A, Chopard B. Theory and applications of an alternative lattice Boltzmann grid refinement algorithm. Physical Review E, 2003, 67(6): 066707 doi: 10.1103/PhysRevE.67.066707
[16] Touil H, Ricot D, Lévêque E. Direct and large-eddy simulation of turbulent flows on composite multi-resolution grids by the lattice Boltzmann method. Journal of Computational Physics, 2014, 256: 220-233 doi: 10.1016/j.jcp.2013.07.037
[17] Rohde M, Kandhai D, Derksen JJ, et al. A generic, mass conservative local grid refinement technique for lattice-Boltzmann schemes. International Journal for Numerical Methods in Fluids, 2006, 51(4): 439-468 doi: 10.1002/fld.1140
[18] Yu DZ, Mei RW, Shyy W. A multi-block lattice Boltzmann method for viscous fluid flows. International Journal for Numerical Methods in Fluids, 2002, 39(2): 99-120 doi: 10.1002/fld.280
[19] Astoul T, Wissocq G, Boussuge JF, et al. Analysis and reduction of spurious noise generated at grid refinement interfaces with the lattice Boltzmann method. Journal of Computational Physics, 2020, 418: 109645 doi: 10.1016/j.jcp.2020.109645
[20] Astoul T, Wissocq G, Boussuge JF, et al. Lattice Boltzmann method for computational aeroacoustics on non-uniform meshes: A direct grid coupling approach. Journal of Computational Physics, 2021, 447: 110667 doi: 10.1016/j.jcp.2021.110667
[21] Peng Y, Shu C, Chew YT, et al. Application of multi-block approach in the immersed boundary–lattice Boltzmann method for viscous fluid flows. Journal of Computational Physics, 2006, 218(2): 460-478 doi: 10.1016/j.jcp.2006.02.017
[22] Chapman S, Cowling TG. The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases. Cambridge: Cambridge University Press, 1990
[23] Huang RZ, Wu HY. Multiblock approach for the passive scalar thermal lattice Boltzmann method. Physical Review E, 2014, 89(4): 043303 doi: 10.1103/PhysRevE.89.043303
[24] d'Humières D. Multiple–relaxation–time lattice Boltzmann models in three dimensions. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 2002, 360(1792): 437-451 doi: 10.1098/rsta.2001.0955
[25] Liou TM, Wang CS. Three-dimensional multidomain lattice Boltzmann grid refinement for passive scalar transport. Physical Review E, 2018, 98(1): 013306 doi: 10.1103/PhysRevE.98.013306
[26] Bhatnagar PL, Gross EP, Krook M. A model for collision processes in gases. i. Small amplitude processes in charged and neutral one-component systems. Physical Review, 1954, 94(3): 511 doi: 10.1103/PhysRev.94.511
[27] 师羊羊. 非平衡流动的高阶格子玻尔兹曼方法. [博士论文]. 哈尔滨: 哈尔滨工业大学, 2021 Shi Yangyang. High order lattice Boltzmann method for non-equilibrium flows. [PhD Thesis]. Harbin: Harbin Institute of Technology, 2021 (in Chinese)
[28] Shan XW, Yuan XF, Chen HD. Kinetic theory representation of hydrodynamics: A way beyond the Navier–Stokes equation. Journal of Fluid Mechanics, 2006, 550: 413-441 doi: 10.1017/S0022112005008153
[29] Guo ZL, Xu K, Wang RJ. Discrete unified gas kinetic scheme for all knudsen number flows: Low-speed isothermal case. Physical Review E, 2013, 88(3): 033305 doi: 10.1103/PhysRevE.88.033305
[30] He XY, Shan XW, Doolen GD. Discrete Boltzmann equation model for nonideal gases. Physical Review E, 1998, 57(1): R13 doi: 10.1103/PhysRevE.57.R13
[31] 李旭晖, 单肖文, 段文洋. 格子玻尔兹曼正则化碰撞模型的理论进展. 空气动力学学报, 2022, 40(3): 1-20 (Li Xuhui, Shan Xiaowen, Duan Wenyang. The theory progress on the regularized lattice Boltzmann collision models. Acta Aerodynamica Sinica, 2022, 40(3): 1-20 (in Chinese) LI Xuhui, Shan Xiaowen, Duan Wenyang. The theory progress on the regularized lattice Boltzmann collision models[J]. Acta Aerodynamica Sinica, 2022, 40(3):1-20 (in Chinese).
[32] Guo ZL, Zheng CG, Shi BC. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Physical Review E, 2002, 65(4): 046308 doi: 10.1103/PhysRevE.65.046308
[33] Krüger T, Varnik F, Raabe D. Second-order convergence of the deviatoric stress tensor in the standard Bhatnagar-Gross-Krook lattice Boltzmann method. Physical Review E, 2010, 82(2): 025701 doi: 10.1103/PhysRevE.82.025701
[34] Lallemand P, Luo LS. Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability. Physical Review E, 2000, 61(6): 6546 doi: 10.1103/PhysRevE.61.6546
[35] Dong ZQ, Wang LP, Peng C, et al. A systematic study of hidden errors in the bounce-back scheme and their various effects in the lattice Boltzmann simulation of viscous flows. Physics of Fluids, 2022, 34(9): 093608 doi: 10.1063/5.0106954
[36] Chai ZH, Shi BC, et al. Multiple-relaxation-time lattice Boltzmann method for the Navier-Stokes and nonlinear convection-diffusion equations: Modeling, analysis, and elements. Physical Review E, 2020, 102(2): 023306
[37] Chen SY, Peng C, Teng YH, et al. Improving lattice Boltzmann simulation of moving particles in a viscous flow using local grid refinement. Computers & Fluids, 2016, 136: 228-246
[38] Shi YY, Wu L, Shan XW. Accuracy of high-order lattice Boltzmann method for non-equilibrium gas flow. Journal of Fluid Mechanics, 2021, 907: A25
[39] Gendre F , Ricot D , Fritz G, et al. Grid refinement for aeroacoustics in the lattice Boltzmann method: A directional splitting approach. Physical Review E, 2017, 96(2): 023311 doi: 10.1103/PhysRevE.96.023311
[40] Chen HD, Filippova O, Hoch J, et al. Grid refinement in lattice Boltzmann methods based on volumetric formulation. Physica A: Statistical Mechanics and its Applications, 2006, 362(1): 158-167 doi: 10.1016/j.physa.2005.09.036
[41] Min HD, Peng C, Guo ZL, et al. An inverse design analysis of mesoscopic implementation of non-uniform forcing in mrt lattice Boltzmann models. Computers & Mathematics with Applications, 2019, 78(4): 1095-1114
[42] Shi BC, Guo ZL. Lattice Boltzmann model for nonlinear convection-diffusion equations. Physical Review E, 2009, 79(1): 016701 doi: 10.1103/PhysRevE.79.016701
[43] Bird RB. Transport phenomena. Appl. Mech. Rev., 2002, 55(1): R1-R4 doi: 10.1115/1.1424298
[44] Flekkøy EG. Lattice bhatnagar-gross-krook models for miscible fluids. Physical Review E, 1993, 47(6): 4247 doi: 10.1103/PhysRevE.47.4247
[45] Guo ZL, Shi BC, Wang NC. Fully lagrangian and lattice Boltzmann methods for the advection-diffusion equation. Journal of Scientific Computing, 1999, 14: 291-300 doi: 10.1023/A:1023273603637
[46] Zhang MX, Zhao WF, Lin P. Lattice Boltzmann method for general convection-diffusion equations: MRT model and boundary schemes. Journal of Computational Physics, 2019, 389: 147-163 doi: 10.1016/j.jcp.2019.03.045
[47] Zheng HW, Shu Chang, Chew YT. A lattice Boltzmann model for multiphase flows with large density ratio. Journal of Computational Physics, 2006, 218(1): 353-371 doi: 10.1016/j.jcp.2006.02.015
[48] Chai ZH, Zhao TS. Lattice Boltzmann model for the convection-diffusion equation. Physical Review E, 2013, 87(6): 063309 doi: 10.1103/PhysRevE.87.063309
[49] Guo ZL, Zheng CG, Shi BC. An extrapolation method for boundary conditions in lattice Boltzmann method. Physics of Fluids, 2002, 14(6): 2007-2010 doi: 10.1063/1.1471914
[50] Zou QS, He XY. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Physics of Fluids, 1997, 9(6): 1591-1598 doi: 10.1063/1.869307
[51] AbdelMigid TA, Saqr KM, Kotb MA , et al. Revisiting the lid-driven cavity flow problem: Review and new steady state benchmarking results using gpu accelerated code. Alexandria Engineering Journal, 2017, 56(1): 123-135 doi: 10.1016/j.aej.2016.09.013
[52] Ladd AJC. Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation. Journal of Fluid Mechanics, 1994, 271: 285-309 doi: 10.1017/S0022112094001771
[53] Nie XB, Shan XW, Chen H. Galilean invariance of lattice Boltzmann models. Europhysics Letters, 2008, 81(3): 34005 doi: 10.1209/0295-5075/81/34005
[54] Wan DD, Wang GC, Chen SY. Numerical investigation of lid-driven deep cavity with local grid refinement of MRT-LBM. Journal of Beijing Institute of Technology, 2019, 28(3): 536-548
[55] Guzik SM, Weisgraber TH, Colella P, et al. Interpolation methods and the accuracy of lattice-Boltzmann mesh refinement. Journal of Computational Physics, 2014, 259: 461-487 doi: 10.1016/j.jcp.2013.11.037
[56] Eitel-Amor G, Meinke M, Schröder W. A lattice-Boltzmann method with hierarchically refined meshes. Computers & Fluids, 2013, 75: 127-139
[57] Lagrava D, Malaspinas O, Latt J, et al. Advances in multi-domain lattice Boltzmann grid refinement. Journal of Computational Physics, 2012, 231(14): 4808-4822 doi: 10.1016/j.jcp.2012.03.015
[58] Grad H. Note on n-dimensional hermite polynomials. Communications on Pure and Applied Mathematics, 1949, 2(4): 325-330 doi: 10.1002/cpa.3160020402
计量
- 文章访问数: 695
- HTML全文浏览量: 181
- PDF下载量: 176
