RESEARCH ON THE LATTICE BOLTZMANN ALGORITHM FOR GRID REFINEMENT BASED ON NON-UNIFORM RECTANGULAR GRID
-
摘要: 传统的格子玻尔兹曼方法(LBM), 特别是基于均匀正方形网格的经典单松弛计算模型(SLBM), 其算法鲁棒性和数值稳定性较差, 限制了LBM的发展和应用. 而网格细化策略可以有效缓解这一窘境, 但是传统LBM中网格细化必然会导致计算效率骤降, 计算设备要求攀高. 为了解决这一问题, 文章基于非均匀矩形网格结构, 结合插值LBM算法的思路, 在保证物面处和流动变化剧烈区域的局部网格细化以及计算精度的前提下, 提出了25点拉格朗日插值LBM算法. 以经典顶盖驱动方腔内流为算例, 开展了包括不同网格分辨率和插值格式的对比分析研究. 验证算例既包括了定常流动的数值模拟, 也涉及了非定常周期性流动的求解. 计算结果表明, 相较于其他插值格式, 拉格朗日插值格式表现优异; 文章局部网格细化工作可以确保物面处及流动变化剧烈区域流动细节的捕捉; 数值模拟算法可以为数值仿真提供可信的计算结果; 同时大幅降低了总网格数量. 因此很大程度上提升了计算效率; 数值模拟方法鲁棒性较好, 适用于包括定常和非定常流动的数值模拟.Abstract: The traditional lattice Boltzmann method (LBM), especially the classic single-relaxation model (SLBM) based on the uniform square grid, has poor robustness and numerical stability, which limits the development and applications of LBM. Grid refinement strategy can effectively alleviate this dilemma, however for the traditional LBM, the grid refinement will inevitably lead to a sudden drop in computational efficiency and a rise in equipment requirements. Therefore, in order to solve this problem, based on the non-uniform rectangular grid, combined with the idea of interpolation LBM, the 25-bit Lagrangian interpolation LBM is proposed on the premise of ensuring the local grid refinement for the surfaces and area with severe flow changes, and the computational accuracy as well. Taking the classic lid-driven cavity flow for instance, a comparative analysis including different grid resolutions and interpolation schemes is performed. The verification includes both the numerical simulations of steady states and unsteady periodic solutions. The results show that the Lagrangian interpolation scheme performs better than other interpolation schemes. In this paper, the local grid refinement is able to ensure the capture of the flow details adjacent to surfaces and in the area of intense flow changes. The numerical algorithm can provide reliable results for numerical simulations. Meanwhile, the total grid number is greatly reduced, as a result the computational efficiency is greatly improved; The numerical simulation method has good robustness and is suitable for numerical simulations for both steady states and unsteady solutions.
-
引 言
格子玻尔兹曼方法[1-3](LBM)在经历几十年的快速发展后已逐渐成为流体力学研究中的重要分支和热点话题. 同时在计算流体力学(CFD)领域中得到广泛应用, 比如不可压缩流动[4]、可压缩流动[5]、多孔介质流[6]、气动声学[7]、多相流[8]、燃烧[9]、渗流[10]等. 由于其经典的碰撞迁移理论完美地契合均匀正方形网格的网格结构, 因此在LBM相关数值模拟中应用非常广泛. 虽然基于均匀正方形网格的单松弛LBM (传统LBM)在数值模拟方面表现不俗, 但由于其鲁棒性和数值稳定性较差, 尤其当松弛因子$ \tau = {{{\text{0}}{\text{.5 + 3}}UL} \mathord{\left/ {\vphantom {{{\text{0}}{\text{.5 + 3}}UL} {(Re\Delta x)}}} \right. } {(Re\Delta x)}} $接近0.5时, 数值稳定性逐渐丧失[1-3], 因此其应用范围基本仅限于低雷诺数Re的数值模拟. 为了应对逐渐增大的雷诺数, 就必须减小网格间距$ \Delta x $, 同时还能满足物面处以及流动变化剧烈区域流动细节的捕捉. 而均匀正方形网格的细化必然导致网格总量的骤然增加, 使得计算效率大幅下降, 且对计算设备的要求也急剧提高. 极大地限制传统LBM的应用与发展.
作者在之前的研究中提出了基于树网格的LBM算法[11-12], 可以在保证计算精度, 满足局部网格细化要求的前提下, 大幅提高传统LBM的计算效率. 此外, 众多学者也开展了传统LBM局部网格细化的算法优化研究[13-27]. 比如有限差分LBM[13-14]、插值补充LBM[15]、多重网格LBM[16-18]、多区块LBM[19-21]、树网格LBM[22-25]以及均匀矩形网格LBM[26-27]等. Mei等[13]以及Guo等[14]借鉴传统的求解N-S方程的手段, 将有限差分方法与LBM结合. 由于控制方程, 即格子玻尔兹曼方程不再通过碰撞迁移演化求解, 转而依靠差分格式逼近. 所以网格结构的使用也变得较为灵活, 随着贴体网格的引入, 在网格整体数量不变的情况下解决了网格物面处的局部加密的问题. 但是贴体网格的使用会带来复杂的坐标变换以及在求解过程中需要构造高精度差分格式, 总体来讲未能有效提升计算效率. 同样地, 基于贴体网格, He等[15]使用的插值补充格式, 依然保留了LBM的精髓, 即分布函数的碰撞迁移过程, 然后通过插值来近似求解碰撞迁移后虚拟点的分布函数, 但计算过程比较复杂且插值格式精度较低. Marvriplis[16]、Patil等[17]和王兴勇等[18]将传统CFD中的多重网格技术应用于LBM的数值模拟中. 其核心思想是通过粗细网格的信息交换, 使得低频误差在粗网格快速降低, 高频误差在细网格得到光滑处理, 从而提升计算效率. 总体来讲, 多区块LBM[19-21]和树网格LBM[22-25]算法基本是殊途同归, 虽然每个学者的实现手段不尽相同, 但是核心思想是一致的. 即对流场内部的不同区域采取局部加密. 这些区域往往是流动变化较为剧烈的区域或者是物面流动区. Zhou[26-27]提出了基于均匀矩形网格的LBM算法, 不同于本文非均匀矩形网格LBM算法, Zhou的计算模型没有插值过程, 用重新构建的矩形离散速度模型以及对应的平衡态分布函数取代传统的正方形离散速度模型(见图1), 可以保证分布函数以矩形方式完成格点之间的碰撞迁移. 虽然该方法可以在一定程度上减少总体网格数量, 但局部网格细化效果不佳且不具备普遍的适用性. 以上国内外的相关研究对本文工作予以很大的启发, 尤其是插值补充LBM算法[15]中虚拟点分布函数的插值求解过程, 对本文基于非均匀矩形网格25点拉格朗日插值LBM算法的构建提供了思路.
1. 25点插值LBM计算模型
本文研究使用最常用的LBGK-D2Q9模型[28]. 其平衡态分布函数$ f_\alpha ^{eq} $的构造如下式所示
$$ f_\alpha ^{eq} = {\omega _\alpha }\rho \left[ {1 + \frac{{{{\boldsymbol{e}}_\alpha } \cdot {\boldsymbol{u}}}}{{c_s^2}} + \frac{{{{\left( {{{\boldsymbol{e}}_\alpha } \cdot {\boldsymbol{u}}} \right)}^2}}}{{2c_s^4}} - \frac{{{{\boldsymbol{u}}^2}}}{{2c_s^2}}} \right] $$ (1) 其中, $ {\omega _\alpha } $是离散速度模型中不同离散方向的权系数; $ \rho $为流体密度; ${\boldsymbol{u}}$为流体速度矢量; ${{\boldsymbol{e}}_\alpha }$是离散速度, 其大小和方向构造如下所示
$$\left.\begin{split} & {{\boldsymbol{e}}_\alpha } = c\left[ \begin{gathered} 0,1,0, - 1,0,1, - 1, - 1,1 \\ 0,0,1,0, - 1,1,1, - 1, - 1 \\ \end{gathered} \right] \\ & {c_s} = \frac{c}{{\sqrt 3 }},\;\;{\omega _\alpha } = \left\{ \begin{gathered} {4 \mathord{\left/ {\vphantom {4 9}} \right. } 9},{\boldsymbol{e}}_\alpha ^2 = 0 \\ {1 \mathord{\left/ {\vphantom {1 9}} \right. } 9},{\boldsymbol{e}}_\alpha ^2 = {c^2} \\ {1 \mathord{\left/ {\vphantom {1 {36}}} \right. } {36}},{\boldsymbol{e}}_\alpha ^2 = 2{c^2} \\ \end{gathered} \right. \\ & \alpha = 0,1,2,\cdots,8 \end{split}\right\} $$ (2) 其中, $ c = {{\Delta x} \mathord{\left/ {\vphantom {{\Delta x} {\Delta t}}} \right. } {\Delta t}} = 1 $为格子速度; $ {c_s} $为格子声速; $ \Delta x $为网格步长; $ \Delta t $为时间步长. 图1介绍了LBGK-D2Q9计算模型的离散速度模型. 详细说明了分布函数$ {f_\alpha } $的碰撞迁移方式.
式(3)定义了LBGK计算模型中碰撞迁移理论的控制方程, 其中${\boldsymbol{ r}}$为空间向量; ${\boldsymbol{\xi}}$为速度向量; $ \tau $ 为松弛时间; $ {f_\alpha } $和$ f_\alpha ^{eq} $分别为分布函数和平衡态的分布函数.
$$ \begin{split} & {f_\alpha }({\boldsymbol{r}} + {{\boldsymbol{e}}_\alpha }\Delta t,t + \Delta t) - {f_\alpha }({\boldsymbol{r}},t) = \\ &\qquad \frac{1}{\tau }\left[ {f_\alpha ^{eq}({\boldsymbol{r}},{\boldsymbol{\xi}} ) - {f_\alpha }({\boldsymbol{r}},{\boldsymbol{\xi}} ,t)} \right]\end{split} $$ (3) 下式定义了流场宏观物理量, 速度${\boldsymbol{u}}$和密度$ \rho $与微观粒子分布函数$ {f_\alpha } $之间的联系. 每一次演化过后, 流场的速度和密度都会通过下式得到更新
$$ \left.\begin{split} & \rho = \sum\limits_\alpha {{f_\alpha }} \\ & \rho {\boldsymbol{u }}= \sum\limits_\alpha {{{\boldsymbol{e}}_\alpha }{f_\alpha }} \end{split} \right\} $$ (4) 传统的LBM算法基于均匀正方形网格, 根据其碰撞迁移理论, 分布函数的迁移过程必须保证在网格格点上. 因为在传统的LBM算法中, 格子速度$ \psi $. 如图2所示, 点A和O为网格格点. 点O处的流体粒子, 其2方向的分布函数经过一个时间步长$ \psi $的演化会迁移至A点, 迁移距离为网格间距$ \omega $, 即$ \mathop {{f_2}}\limits^A = \mathop {{f_2}}\limits^O $. 其他方向的分布函数也遵循同样的规则迁移至O点周围的各点. 由此可见均匀正方形网格天然匹配LBM的演化机理. 通常情况下, 为了得到较为理想的数值模拟结果, 尤其是捕捉物面处以及流动剧烈变化区的流动细节, 往往需要极为细腻的网格. 这就使得均匀正方形网格的网格总量急剧增加, 极大地降低了传统LBM的计算效率和鲁棒性, 限制了LBM的应用. 因此, 本文采用了非均匀矩形网格用于流场局部加密既, 保证了物面处和流动变化剧烈区的流场细节捕捉, 同时相较于均匀正方形网格, 在保证计算精度的前提下又能提升LBM的计算效率.
由于本文的计算网格是非均匀的矩形网格, 每一层网格的网格间距都不相同, 传统LBM中基于格点和格点之间迁移条件无法满足.
如图3所示, 点A和O均为网格格点. 由于点A所在的网格层和点O所在的网格层之间的网格间距不同, $ \Delta {x_{OA}} \lt \Delta {x_{OC}} $. 所以当经过一个时间步长$ \Delta t $的演化后, 针对O点, 只有2方向的分布函数$ {f_2} $可以正常迁移到同为网格格点的A点($ \mathop {{f_2}}\limits^A = \mathop {{f_2}}\limits^O $). 其他方向的分布函数均无法完成正常迁移, 比如, 7方向的分布函数只能迁移至B'点, 而此点为虚拟点, 且我们并不知道此点的信息. 因此本文数值模拟, 除了LBM原有的碰撞迁移计算还需要额外构造一个插值格式近似求解B'点的信息. 在插值LBM算法[15]的基础上, 本文提出了基于25点拉格朗日插值格式的非均匀矩形网格LBM算法.
图4, 以$ {A_{33}} $点5方向的分布函数$f_5$为例, 介绍了本文25点拉格朗日插值格式的构造过程. 由图可见点${A_{ij}}\;(i=1,2,\cdots,5;j=1,2,\cdots,5)$均为计算域的实际网格点, 而点$ {a_{ij}} $($i=1,2,\cdots,5;j=1,2,\cdots,5$)为沿5方向迁移后的虚拟点, 并且包括点$ {A_{33}} $以及其周围的24个点都参与了插值格式的构建. 网格沿横(${{\xi }_{1}},{{\xi }_{2}},{{\xi }_{3}}$和${{\xi }_{4}}$)、纵(${{\eta }_{1}},{{\eta }_{2}},{{\eta }_{3}}$和${{\eta }_{4}}$)向均有4层网格且网格间距均不相同. $\varDelta = \Delta t \cdot c$为流体粒子单位时间步长$ \Delta t $内移动的距离($ c = 1.0 $), 同时也是最细网格间距($\varDelta = \Delta {x_{{\rm{finest}}}}$). 本文沿横、纵向靠近物面的第一层最细网格为$\Delta {{x}_{{\rm{finest}}}} = 1/1024$, 见图5(b)中计算域四角黑色区域. 以5方向的分布函数为例, 下式给出了本文25点拉格朗日插值格式
![]()
图 4 非均匀矩形网格的迁移演示($ \alpha = 5 $)Figure 4. Streaming on discrete direction ($ \alpha = 5 $)$$ \mathop {{f_5}}\limits^{{A_{33}}} = \sum\limits_{j = 1}^5 {\sum\limits_{i = 1}^5 {{\gamma _j}{\varsigma _i}\mathop {{f_5}}\limits^{{{{a}}_{ij}}} } } $$ (5) 其中, ${{\gamma }_{j}}$和${\varsigma _i}$是本文25点拉格朗日插值格式的系数. 而$ \mathop {{f_5}}\limits^{{a_{ij}}} $是虚拟点$ {a_{ij}} $沿5方向从点$ {A_{ij}} $迁移过来的分布函数($ \mathop {{f_5}}\limits^{{a_{ij}}} = \mathop {{f_5}}\limits^{{A_{ij}}} $).
2. 数值模拟
2.1 计算域构建及边界条件
为了验证本文基于25点拉格朗日插值的非均匀矩形网格LBM算法, 选取经典的顶盖驱动方腔内流为基准算例($Ma = 0.173\;2$且$Re < {\text{10 0}}00$). 本文数值模拟最细网格步长($\Delta {{x}_{{\rm{finest}}}} = 1/1024$)的选取借鉴了作者之前研究工作[29-32]中针对腔体内流的网格独立性验证. 基于此网格尺度, 发现当$Re < {\text{1}}0\;000$时能够保证$ {y^ + } \lt {\text{0}}{\text{.1}} $, 数值模拟结果可信且能够捕捉边界层附近的流场信息. 如图5所示, 展示了顶盖驱动方腔内流的网格结构, 图5(a)和图5(b)分别对应均匀正方形网格和非均匀矩形网格. 流动信息采集点位于腔体中心$ P(x = {L \mathord{\left/ {\vphantom {L 2}} \right. } 2},y = {L \mathord{\left/ {\vphantom {L 2}} \right. } 2}) $处.
为了便于研究, 本文非均匀矩形网格的建立借鉴了等比数列的控制函数, 将计算域沿横向和纵向划分为4个对称的子区间. 每个子区间网格间距由下式控制
$$\left. \begin{split} & \Delta {x_i} = \Delta {x_{{\rm{finest}}}}{(1 + \sigma )^{i - 1}},{\text{ }}i = 0,1,2,\cdots,{N \mathord{\left/ {\vphantom {N 2}} \right. } 2} \\ & \Delta {y_j} = \Delta {y_{{\rm{finest}}}}{(1 + \sigma )^{j - 1}},{\text{ }}j = 0,1,2,\cdots,{N \mathord{\left/ {\vphantom {N 2}} \right. } 2} \\ & \Delta {x_i} = \Delta {x_{{\rm{finest}}}}{(1 + \sigma )^{ - (i - 1)}},{\text{ }}i = {N \mathord{\left/ {\vphantom {N 2}} \right. } 2}{\text{ + 1}},\cdots,N \\ & \Delta {y_j} = \Delta {y_{{\rm{finest}}}}{(1 + \sigma )^{ - (j - 1)}},{\text{ }}j = {N \mathord{\left/ {\vphantom {N {2{\text{ + 1}}}}} \right. } {2{\text{ + 1}}}},\cdots,N \end{split}\right\} $$ (6) 其中, N为网格点个数, $ \sigma $是等比控制函数的公比, 在本文中当$ N = 500 $时, $ \sigma $为0.005 213 64. $\Delta {x_{{\rm{finest}}}} = \Delta {y_{{\rm{finest}}}}$定义了贴近物面处第一层网格的网格间距(最细网格).
2.2 边界条件处理
本文涉及的边界均为平直边界(方腔4个边). 其中顶盖为驱动边界, 其余三边均为物面边界. 由于非平衡态外推格式[33]易于代码编写且非常适用于平直边界点的整体处理. 本文延续之前研究中已成熟应用的非平衡态外推格式, 将边界点上的分布函数分为平衡态($ f_\alpha ^{eq} $)和非平衡态($ f_\alpha ^{neq} $)两部分. 其中平衡态部分由平衡态分布函数的定义近似获得, 而非平衡态部分则用非平衡态外推法求解.
如图6所示点A, B和C为物面边界点, D点为流场点. 根据LBM的演化原理可知在每次演化之前需要知道每个点的分布函数, 对于B点, 其分布函数可分为平衡态和非平衡态两部分
$$ \mathop {{f_\alpha }}\limits^B = \mathop {f_\alpha ^{eq}}\limits^B + \mathop {f_\alpha ^{neq}}\limits^B $$ (7) 以物面边界点B为例, 该点的平衡态分布函数$ \mathop {f_\alpha ^{eq}}\limits^B $可用该点的宏观物理量来构造, 如果该点的宏观物理量未知, 则由D点的相应值代替. 而非平衡态分布函数则由D点的非平衡态分布函数来近似代替
$$ \mathop {f_\alpha ^{neq}}\limits^B = \mathop {{f_\alpha }}\limits^D - \mathop {f_\alpha ^{eq}}\limits^D $$ (8) 从而得到B点的分布函数$ \mathop {{f_\alpha }}\limits^B $
$$ \mathop {{f_\alpha }}\limits^B = \mathop {f_\alpha ^{eq}}\limits^B + \mathop {{f_\alpha }}\limits^D - \mathop {f_\alpha ^{eq}}\limits^D $$ (9) 3. 计算结果与分析讨论
3.1 定常流动数值模拟
本文以顶盖驱动方腔内流为例, 针对本文RLBM-L25算法在定常($ Re = 1000 $和$ Re = {\text{5}}000 $)和非定常状态($ Re = {\text{88}}00 $和$ Re = {\text{95}}00 $)下开展数值模拟分析验证研究.
表1介绍了当$ Re = 1000 $时本文的计算结果和其他文献[34-35]的数据对比. 第1列P代表流场中涡漩的位置, LB表示左下侧次级涡; RB表示右下侧次级涡; PC代表中心主涡(见图7). 对应每个涡结构给出了4个流动参数, 分别为涡心横坐标(X)、纵坐标(Y)、流函数($ \psi $)和涡量($ \omega $). 由表1可知, 本文提出的RLBM-L25算法其结果与其他文献的数据吻合较好, 可以输出可信的数值模拟结果.
表 1 不同计算结果对比 (${\boldsymbol{Re}} = {\boldsymbol{1000}}$)Table 1. Different results comparison ($ Re = 1000 $)P Ghia et al.[34] Erturk et al.[35] This work LB X = 0.085 9 X = 0.083 3 X = 0.082 4 Y = 0.078 1 Y = 0.078 3 Y = 0.07698 $ \psi $ = 0.000231 $ \psi $ = 0.000233 $ \psi $ = 0.000228 $ \omega $ = 0.36175 $ \omega $ = 0.35427 $ \omega $ = 0.346 2 RB X = 0.859 4 X = 0.863 3 X = 0.863 5 Y = 0.109 4 Y = 0.111 7 Y = 0.111 4 $ \psi $ = 0.001751 $ \psi $ = 0.00173 $ \psi $ = 0.001716 $ \omega $ = 1.154 7 $ \omega $ = 1.118 2 $ \omega $ = 1.092 1 CP X = 0.531 3 X = 0.530 0 X = 0.527 9 Y = 0.562 5 Y = 0.565 0 Y = 0.562 5 $ \psi $ = −0.11793 $ \psi $ = −0.118 9 $ \psi $ = −0.11813 $ \omega $ = −2.04968 $ \omega $ = −2.067 8 $ \omega $ = −2.047 2 类似表1, 表2列出了当$ Re = {\text{5}}000 $时的数据对比. 数据结构和内容与表1一致. 再次证明RLBM-L25算法对于定常流动数值模拟表现良好.
表 2 不同计算结果对比 (${\boldsymbol{Re}} = {\boldsymbol{5000}}$)Table 2. Different results comparison ($ Re = {\text{5}}000 $)P Ghia et al.[34] Erturk et al.[35] This work LB X = 0.070 3 X = 0.073 3 X = 0.073 0 Y = 0.136 7 Y = 0.136 7 Y = 0.135 3 $ \psi $ = 0.001361 $ \psi $ = 0.001376 $ \psi $ = 0.001383 $ \omega $ = 1.530 6 $ \omega $ = 1.514 3 $ \omega $ = 1.521 RB X = 0.808 6 X = 0.805 0 X = 0.803 2 Y = 0.074 2 Y = 0.073 3 Y = 0.072 5 $ \psi $ = 0.003084 $ \psi $ = 0.003073 $ \psi $ = 0.003089 $ \omega $ = 2.663 5 $ \omega $ = 2.739 $ \omega $ = 2.743 CP X = 0.511 7 X = 0.515 0 X = 0.512 0 Y = 0.535 2 Y = 0.535 0 Y = 0.532 1 $ \psi $ = −0.119 0 $ \psi $ = −0.122 2 $ \psi $ = −0.119 8 $ \omega $ = −1.860 2 $ \omega $ = −1.940 5 $ \omega $ = −1.892 如图7(a)和图7(b)所示, 展示了定常状态下的流场拓扑结构. 可以明显观察到经典方腔驱动内流的典型流动特征, 与其他文献[34-39]的数据完全一致. 中心主涡覆盖大部分计算域从而主导整个流场, 随着雷诺数的增加, 次级涡逐渐增大并伴随边角处出现细碎涡的现象. 图7(c)展示了$ Re = {\text{1}}000 $时流场的局部放大图, 左侧是基于均匀正方形网格的传统LBM算法(SLBM)的计算结果, 右侧是本文RLBM-L25的计算结果. 通过对比可以看出, 本文基于非均匀矩形网格的网格加密方法可以保证物面处和流动突变区的流动细节捕捉. 图7(d)和图7(e)分别展示了$ Re = {\text{1}}000 $时基于均匀正方形网格的经典单松弛LBM(SLBM)且网格分辨率为1024和基于非均匀矩形网格插值LBM(RLBM-500)且网格分辨率为500时的涡量图对比, 如图所示, 涡量图吻合得较好.
图8描述了沿顶盖分布的剪切力$ {f_s} $, 总体来讲, 本文RLBM-L25算法网格分辨率分别为500和600且雷诺数等于1000时的计算结果与传统LBM算法网格分辨率为1024时的结果较为吻合, 但是在顶盖左右两侧附近仍存在一定量级的误差且不会因为插值格式的改变而消除, 因此作者认为这种情况的出现首先跟顶盖上网格分布数量有关, 其次跟网格局部细化的控制方程也有关.
图9描述了$ Re = {\text{5}}000 $时水平速度分量$ {u_x} $的均方根误差${e_{{\rm{RMS}}}} = \sqrt {\dfrac{1}{n}\displaystyle\sum\limits_{i = 1}^n {{{(u_x^i - \bar U)}^2}} }$随时间的变化曲线. 其中$ u_x^i $为第i时步的计算结果; $ \bar U $为最终收敛解. 红色曲线代表本文算法(RLBM-L25-500)的计算结果. 黑色曲线表示传统算法(SLBM-1024)的计算结果. 如图所示, 本文算法的收敛速度快于传统算法, 从而在满足局部网格细化且保证计算精度的前提下提高了LBM的计算效率. 研究发现, 当计算状态和设备相同的情况下, 根据不同的网格细化策略, 可以将计算效率提升1.7 ~ 5.2倍. 本文RLBM-L25算法的收敛速度比传统算法(SLBM)快$ t = N \times \Delta t = 2.6{\text{36}} $倍. 其中$ t $为计算总耗时比; $ N = 1.56 $为收敛所需演化步数比; $ \Delta t = 1.69 $为单位演化步所需的物理计算时间比.
![]()
图 9 速度均方根误差曲线Figure 9. History of root mean square error${e_{{\rm{RMS}}}}$of velocity$ {u_x} $图10展示了扰动衰减率$\varepsilon = \lg \left| {{u_x} - \bar U} \right|$随时间的变化曲线, 由图可见红色曲线线性段明显陡于黑色曲线, 说明相较于SLBM算法, 本文RLBM-L25算法的收敛速度更快, 意味着计算效率更高.
![]()
图 10 扰动衰减率曲线Figure 10. History of perturbation decay rate$ \varepsilon $of velocity$ {u_x} $为了进一步讨论分析本文算法的表现, 还开展了基于不同网格分辨率的研究. 以传统算法SLBM且网格分辨率为1024时的计算结果作为对比.
由表3可见, 第1列给出了不同的网格分辨率(240 ~ 620). 第2列介绍了当$ Re = {\text{1}}000 $时, 信息采集点水平速度分量$ {u_x} $的收敛值. 第3列展示了每个网格分辨率相对于SLBM-1024算法计算结果的相对误差$ \Delta {u_x}(\% ) $. 由表可见, 随着网格分辨率的增加, 相对误差逐渐减小, 这与均匀正方形网格数值模拟所表现出来的特性一致. 当网格分辨率从500增加至620时, 相对误差从1.7%减小至0.99%. 但是由于网格数量的增加会降低计算效率, 同时1.7%的相对误差小于2.0%, 是本文可接受的误差范围. 综合考虑, 本文所有的数值模拟都采用500的网格分辨率.
表 3 不同网格数量计算结果对比Table 3. Comparison of different resolutionsResolution $ {u_x} $ $\Delta {u_x}/\%$ RLBM-240 −0.00574936 6.72 RLBM-300 −0.00587757 4.85 RLBM-400 −0.00599459 2.74 RLBM-500 −0.0060579 1.71 RLBM-620 −0.00610233 0.99 SLBM-1024 −0.00616359 0.0 此外, 本文还开展了不同插值格式之间的对比研究. 同样地我们以传统算法SLBM且网格分辨率为1024时的计算结果作为对比. 结合表4, 可以看到第1列列举了不同的插值格式, 分别为拉格朗日9点插值(RLBM-L9)、牛顿9点插值(RLBM-N9)、埃尔米特9点插值(RLBM-H9)以及拉格朗日25点插值(RLBM-L25). 第2列介绍了当$ Re = {\text{1}}000 $时, 信息采集点水平速度分量$ {u_x} $的收敛值. 第3列展示了每个插值格式相对于SLBM-1024计算结果的相对误差$ \Delta {u_x}(\% ) $. 由表可见, 本文RLBM-L25算法相对于其他插值算法表现优异.
表 4 不同插值格式计算结果对比Table 4. Comparison of different interpolationsInterpolation $ {u_x} $ $\Delta {u_x}/\%$ RLBM-L9 −0.00598082 2.96 RLBM-N9 −0.00598082 2.96 RLBM-H9 −0.00597826 3.0 RLBM-L25 −0.0060579 1.71 SLBM −0.00616359 0.0 4.2 非定常流动分数值模拟
为了进一步探讨本文算法(RLBM-L25)对于非定常流动的表现, 本文开展了非定常周期性流动的数值模拟. 选取两个非定常周期性雷诺数8800和9500作为数值模拟计算条件.
表5展示了不同网格分辨率(RLBM-L25-500和RLBM-L25-620)雷诺数为9500时计算结果, 对SLBM-1024的计算结果. 第1列数据(Data)分别列举了周期性速度曲线的最大值($ {u_{x\max }} $)、最小值($ {u_{x\min }} $)、平均值(${u_{x{\rm{mean}}}}$)以及平均值相对SLBM-1024计算结果的相对误差$\Delta {u_{x{\rm{mean}}}}$. 通过表5我们发现, 对于非定常周期性流动, 随着网格分辨率的增大, 相对误差逐渐减小. 相对于网格分辨率620, 500的分辨率产生较大的相对误差, 但仍然在误差允许范围内($\Delta {u_{x{\rm{mean}}}} \lt 5\%$). 作者认为, 这个网格分辨率可以保证非定常周期性流动的数值模拟精度.
表 5 不同网格分辨率计算结果对比Table 5. Comparison of different resolutionsData RLBM-500 RLBM-620 SLBM-1024 $ {u_{x\max }} $ −0.0024902 −0.00251714 −0.0025769 $ {u_{x\min }} $ −0.0025351 −0.00256417 −0.002624 ${u_{x{\rm{mean}}} }$ −0.0025127 −0.00254066 −0.0026005 $\Delta {u_{x{\rm{mean}}} }$/% 3.37 2.3 0.0 如图11所示, 介绍了流场在一个完整周期($T = {1 \mathord{\left/ {\vphantom {1 f}} \right. } f},f = 0.432\;6$)内某一时刻的瞬时流线图. 其表现出的流场典型特征与作者之前针对方腔内流的研究[41-42]完全一致.
图12介绍了雷诺数为9500时, 信息采集点基于不同插值格式的计算结果对比. 主画面展示了速度频谱图对比, 可见不同的插值格式对于周期性流动的振荡频率(居中内嵌图)的预测基本一致, 且本文RLBM-L25算法表现优异. 从右侧内嵌图可以看出, 在预测周期性流动的速度振幅和平均值时, 与SLBM-1024的计算结果存在差异. 相较于其他插值格式, RLBM-L25算法的表现最好, 即便如此, 与SLBM-1024的计算结果的相对误差为3.5%左右. 作者认为导致这一误差的可能原因有两个. 其一, 与本文网格细化策略有关, 比如网格分辨率, 网格间距控制函数等. 此外插值格式也是一个需要考虑的因素, 作者认为更精确的插值格式能很大程度上消除或降低这一误差.
5. 结 论
本文提出了一种基于25点拉格朗日插值的非均匀矩形网格LBM算法. 针对本文算法以方腔顶盖驱动内流为例开展了定常和非定常流动的数值模拟验证研究, 得出主要结论如下.
(1) 本文算法对于定常流动数值模拟表现优异, 网格细化策略可以保证物面处和流动突变区的流动细节的捕捉.
(2) 本文算法对于非定常周期性流动数值模拟表现良好, 同样可以保证物面处和流动突变区的流动细节的捕捉. 相对基于均匀正方形网格($ 1024 \;\times\; 1024 $)传统LBM的计算结果, 虽然存在误差, 但总体控制在可接受的范围内.
(3) 相较于传统的基于均匀正方形网格的LBM算法, 本文算法在保证计算精度的前提下可以提升计算效率2 ~ 5倍.
(4) 针对本文非均匀矩形网格, 当网格分辨率逐渐增加时, 基于理论解的相对误差逐渐减小, 与均匀正方形网格LBM算法的特性一致.
(5) 本文拉格朗日25点插值格式, 相较于格式拉格朗日9点插值、牛顿9点插值和埃尔米特9点插值, 表现更好.
(6) 为了消除或降低误差, 可以改良网格细化策略, 如提高网格分辨率和优化网格间距控制函数等. 此外, 采用更为精确的插值格式.
-
![]()
图 4 非均匀矩形网格的迁移演示($ \alpha = 5 $)
Figure 4. Streaming on discrete direction ($ \alpha = 5 $)
![]()
图 9 速度均方根误差曲线
Figure 9. History of root mean square error${e_{{\rm{RMS}}}}$of velocity$ {u_x} $
![]()
图 10 扰动衰减率曲线
Figure 10. History of perturbation decay rate$ \varepsilon $of velocity$ {u_x} $
表 1 不同计算结果对比 (${\boldsymbol{Re}} = {\boldsymbol{1000}}$)
Table 1 Different results comparison ($ Re = 1000 $)
P Ghia et al.[34] Erturk et al.[35] This work LB X = 0.085 9 X = 0.083 3 X = 0.082 4 Y = 0.078 1 Y = 0.078 3 Y = 0.07698 $ \psi $ = 0.000231 $ \psi $ = 0.000233 $ \psi $ = 0.000228 $ \omega $ = 0.36175 $ \omega $ = 0.35427 $ \omega $ = 0.346 2 RB X = 0.859 4 X = 0.863 3 X = 0.863 5 Y = 0.109 4 Y = 0.111 7 Y = 0.111 4 $ \psi $ = 0.001751 $ \psi $ = 0.00173 $ \psi $ = 0.001716 $ \omega $ = 1.154 7 $ \omega $ = 1.118 2 $ \omega $ = 1.092 1 CP X = 0.531 3 X = 0.530 0 X = 0.527 9 Y = 0.562 5 Y = 0.565 0 Y = 0.562 5 $ \psi $ = −0.11793 $ \psi $ = −0.118 9 $ \psi $ = −0.11813 $ \omega $ = −2.04968 $ \omega $ = −2.067 8 $ \omega $ = −2.047 2 
下载: 导出CSV
表 2 不同计算结果对比 (${\boldsymbol{Re}} = {\boldsymbol{5000}}$)
Table 2 Different results comparison ($ Re = {\text{5}}000 $)
P Ghia et al.[34] Erturk et al.[35] This work LB X = 0.070 3 X = 0.073 3 X = 0.073 0 Y = 0.136 7 Y = 0.136 7 Y = 0.135 3 $ \psi $ = 0.001361 $ \psi $ = 0.001376 $ \psi $ = 0.001383 $ \omega $ = 1.530 6 $ \omega $ = 1.514 3 $ \omega $ = 1.521 RB X = 0.808 6 X = 0.805 0 X = 0.803 2 Y = 0.074 2 Y = 0.073 3 Y = 0.072 5 $ \psi $ = 0.003084 $ \psi $ = 0.003073 $ \psi $ = 0.003089 $ \omega $ = 2.663 5 $ \omega $ = 2.739 $ \omega $ = 2.743 CP X = 0.511 7 X = 0.515 0 X = 0.512 0 Y = 0.535 2 Y = 0.535 0 Y = 0.532 1 $ \psi $ = −0.119 0 $ \psi $ = −0.122 2 $ \psi $ = −0.119 8 $ \omega $ = −1.860 2 $ \omega $ = −1.940 5 $ \omega $ = −1.892
下载: 导出CSV
表 3 不同网格数量计算结果对比
Table 3 Comparison of different resolutions
Resolution $ {u_x} $ $\Delta {u_x}/\%$ RLBM-240 −0.00574936 6.72 RLBM-300 −0.00587757 4.85 RLBM-400 −0.00599459 2.74 RLBM-500 −0.0060579 1.71 RLBM-620 −0.00610233 0.99 SLBM-1024 −0.00616359 0.0
下载: 导出CSV
表 4 不同插值格式计算结果对比
Table 4 Comparison of different interpolations
Interpolation $ {u_x} $ $\Delta {u_x}/\%$ RLBM-L9 −0.00598082 2.96 RLBM-N9 −0.00598082 2.96 RLBM-H9 −0.00597826 3.0 RLBM-L25 −0.0060579 1.71 SLBM −0.00616359 0.0
下载: 导出CSV
表 5 不同网格分辨率计算结果对比
Table 5 Comparison of different resolutions
Data RLBM-500 RLBM-620 SLBM-1024 $ {u_{x\max }} $ −0.0024902 −0.00251714 −0.0025769 $ {u_{x\min }} $ −0.0025351 −0.00256417 −0.002624 ${u_{x{\rm{mean}}} }$ −0.0025127 −0.00254066 −0.0026005 $\Delta {u_{x{\rm{mean}}} }$/% 3.37 2.3 0.0
下载: 导出CSV
-
[1] Chen SY, Doolen G. Lattice Boltzmann method for fluid flows. Annual Review of Fluid Mechanics, 1998, 143(2): 426-448
[2] Benzi R, Succi S, Vergassola M, The lattice Boltzmann equation: theory and applications. Physics Reports, 1992, 222(3): 145-197
[3] 何雅玲, 王勇, 李庆. 格子Boltzmann方法的理论及应用. 北京: 科学出版社, 2009: 145-197 He Yaling, Wang Yong, Li Qing. Lattice Boltzmann Method: Theory and Applications. Beijing: Science Press, 2009: 145-197 (in Chinese)
[4] Silva G, Semiao V. First- and second-order forcing expansions in a lattice Boltzmann method reproducing isothermal hydrodynamics in artificial compressibility form. Journal of Fluid Mechanics, 2012, 698: 282-303
[5] Yang LM, Shu C, Wu J. Development and comparative studies of three non-free parameter lattice Boltzmann models for simulation of compressible flows. Advances in Applied Mathematics and Mechanics, 2012, 4(4): 454-472
[6] Jin Y, Uth MF, Kuznetsov AV, et al. Numerical investigation of the possibility of macroscopic turbulence in porous media: a direct numerical simulation study. Journal of Fluid Mechanics, 2015, 766: 76-103
[7] Li XM, Leung RCK, So RMC. One-step aeroacoustics simulation using lattice Boltzmann method. AIAA Journal, 2006, 44(1): 78-79
[8] Zhang CY, Cheng P, Minkowycz WJ. Lattice Boltzmann simulation of forced condensation flow on a horizontal cold surface in the presence of a non-condensable gas. International Journal of Heat and Mass Transfer, 2017, 115: 500-512
[9] Chen S, Liu ZH, He Z, et al. A new numerical approach for fire simulation. International Journal of Modern Physics C, 2007, 118(2): 187-202
[10] Haddadi H, Morris JF. Microstructure and rheology of finite inertia neutrally buoyant suspensions. Journal of Fluid Mechanics, 2014, 749: 431-459
[11] An B, Bergadà JM, Mellibovsky F, et al. New Applications of numerical simulation based on lattice Boltzmann method at high Reynolds numbers. Computers & Mathematics with Applications, 2020, 79(6): 1718-1741
[12] 安博, 桑为民. 基于不同网格结构的LBM研究. 力学学报, 2013, 45(5): 699-706 (An Bo, Sang Weimin. The numerical study of lattice Boltzmann method based on different grid structure. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(5): 699-706 (in Chinese) An Bo, Sang Weimin. The numerical study of lattice Boltzmann method based on different grid structure. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(5): 699-706(in Chinese))
[13] Mei RW, Shyy W. On the finite difference-based lattice Boltzmann method in curvilinear coordinates. Journal of Computational Physics, 1998, 143(2): 426-448
[14] Guo ZL, Zhao TS. Explicit finite-difference lattice Boltzmann method for curvilinear coordinates. Physical Review E, 2003, 67(6): 066709
[15] He XY, Doolen G. Lattice Boltzmann method on curvilinear coordinates system: Flow around a circular cylinder. Journal of Computational Physics, 1997, 134(2): 306-315
[16] Marvriplis DJ. Multigrid solution of the steady state lattice Boltzmann equation. Computers & Fluids, 2006, 35: 793-804
[17] Patil DV, Premnath KN, Banerjee S. Multigrid lattice Boltzmann method for accelerated solution of elliptic equations. Journal of Computational Physics, 2014, 265: 172-194
[18] 王兴勇, 索丽生, 程永光等. 双重网格lattice Boltzmann方法. 海河大学学报, 2003, 31(1): 5-10 (Wang Xingyong, Suo Lisheng, Cheng Yongguang, et al. Lattice Boltzmann method with double meshes. Journal of Hohai University, 2003, 31(1): 5-10 (in Chinese) Wang Xingyong, Suo Lisheng, Cheng Yongguang, et al. . Lattice Boltzmann method with double meshes. Journal of Hohai University, 2003, 31(1): 5-10(in Chinese))
[19] Zhang Y, Xie JH, Li XY, et al. A multi-block adaptive solving technique based on lattice Boltzmann method. Modern Physics Letters B, 2018, 32(12-13): 1840052
[20] Lagrava D, Malaspinas O, Latt J, et al. Advances in multi-domain lattice Boltzmann grid refinement. Journal of Computational Physics, 2012, 231: 4808-4822
[21] Lin CL, Lai TG. Lattice Boltzmann method on composite grids. Physical Review E, 2000, 62(2): 2219-2225
[22] Eitel-Amor G, Meinke M, Schrö der W. A lattice-Boltzmann method with hierarchically refined meshes. Computers & Fluids, 2013, 75: 127-139
[23] Lu ZY, Liao Y, Qian DY, et al. Large eddy simulations of a stirred tank using the lattice Boltzmann method on a nonuniform grid. Journal of Computational Physics, 2002, 181: 675-704
[24] Liu B, Khalili A. Acceleration of steady-state lattice Boltzmann simulations for exterior flows. Computers & Fluids, 2013, 75: 127-139
[25] Valero-Lara P, Jansson J. A non-uniform staggered Cartesian grid approach for lattice-Boltzmann method. Procedia Computer Science, 2015, 51: 296-305
[26] Zhou JG. A lattice Boltzmann method for solute transport. International Journal for Numerical Methods in Fluids, 2009, 61: 848-863
[27] Zhou JG. A rectangular lattice Boltzmann method for groundwater flows. Procedia Computer Science, 2007, 21(9): 531-542
[28] Qian YH, d’Humieres D, Lallemand P. Lattice BGK models for Navier-Stokes equation. Europhysics Letters, 1992, 17(6): 478-484
[29] An B, Bergadà JM, Mellibovsky F. The lid driven right-angled isosceles triangular cavity flow. Journal of Fluid Mechanics, 2019, 875: 476-519
[30] An B, Mellibovsky F, Bergadà JM, et al. Towards a better understanding of wall-driven square cavity flows using lattice Boltzmann method. Applied Mathematical Modelling, 2020, 82: 469-486
[31] 安博, 孟欣雨, 桑为民. 镜像对称顶盖驱动方腔内流过渡流临界特性研究. 力学学报, 2022, 54(9): 2409-2418 (An Bo, Meng Xinyu, Sang Weimin. On the transitional characteristics of mirror symmetry lid-driven cavity flow. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(9): 2409-2418 (in Chinese) An Bo, Meng Xinyu, Sang Weimin. On the transitional characteristics of mirror symmetry lid-driven cavity flow. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(9): 2409-2418(in Chinese))
[32] An B, Guo SP, Bergadà JM. Lid driven triangular and trapezoidal cavity flow:Vortical structures for steady solutions and Hopf bifurcations. Applied Sciences, 2023, 13(2): 888
[33] Guo ZL, Zheng CG, Shi BC. An extrapolation method for method boundary conditions in lattice Boltzmann method. Physics of Fluids, 2002, 14(6): 2007-2010
[34] Ghia U, Ghia KN, Shin CT. High-Resolutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, 1982, 48(3): 387-411
[35] Erturk E, Gökçöl C. Fourth-order compact formulation of Navier-Stokes equations and driven cavity flow at high Reynolds numbers. International Journal for Numerical Methods in Fluids, 2005, 50(4): 421-436
[36] Hou S, Zhou Q, Chen S, et al. Simulation of cavity flow by the lattice Boltzmann method. Journal of Computational Physics, 1995, 118(2): 329-347
[37] Vanka S. Block implicit multigrid solution of Navier-Stokes equations in primitive variables. Journal of Computational Physics, 1986, 65(1): 138-158 doi: 10.1016/0021-9991(86)90008-2
[38] Das MK, Rajesh Kanna P. Application of an ADI scheme for steady and periodic solutions in a lid-driven cavity problem. Journal of Numerical Methods for Heat & Fluid Flow, 2007, 17(8): 799-822
[39] Shi X, Huang XW, Zheng Y, et al. A hybrid algorithm of lattice Boltzmann method and finite difference-based lattice Boltzmann method for viscous flows. International Journal for Numerical Methods in Fluids, 2017, 85(11): 641-661
计量
- 文章访问数:
- HTML全文浏览量:
- PDF下载量:
