陈颂英, 汪超, 曲延鹏, 王润堃
山东大学 高效洁净机械制造教育部重点实验室, 山东 济南 250061
收稿日期：2017-05-08
基金项目：山东省科技发展计划项目(2016GGX104018)。
作者简介：陈颂英(1966-)，男，山东莱芜人，山东大学教授，博士生导师。

摘要：使用多参数弛豫模型的格子玻尔兹曼方法及D2Q9模型, 对方腔流顶盖下方左右奇异角落处应用局部网格加密模拟方腔流动特性.在粗细网格界面上, 建立了分布函数的转换公式, 采用3次样条插值计算加密界面上插值点的参量, 根据Chapman-Enskog分析, 给出了应力计算公式.数值计算结果表明:对雷诺数1000的方腔流动, 沿腔中心线的速度分布与经典文献结果对比效果良好, 压力、涡量轮廓图的噪声明显降低, 应力振荡明显减少, 证实了所提局部加密方法的有效性.
关键词：格子Boltzmann多参数弛豫模型局部网格加密方腔D2Q9模型
Local Grid Refinement Implementation of Lid Driven Cavity Flow Based on MRT-LBM
CHEN Song-ying, WANG Chao, QU Yan-peng, WANG Run-kun
Key Laboratory of High-Efficiency and Clean Mechanical Manufacture, Shandong University, Jinan 250061, China
Corresponding author: CHEN Song-ying, E-mail: chensy66@sdu.edu.cn
Abstract: Employing the D2Q9 model, grids at the double singular corners were locally refined to investigate the flow peculiarity of the lid driven cavity based on the multi-relaxation time of lattice Boltzmann method (MRT-LBM). The transferring formulae of the distribution functions were established between the coarse and fine grids. Parameters were interpolated with cubic spline implementation on the refined layers, and the stress expressions were given on the basis of the Chapman-Enskog. For the cavity flow with the Reynolds number 1000, the simulation results showed that the velocity distribution along the cavity mid plane agrees with that of the classical benchmark very well, the noises on the pressure and vortex contours are significantly reduced, and the stress fluctuation obviously decreases. The computational results showed validation of the present method.
Key Words: lattice Boltzmannmulti-relaxation time(MRT)local grid refinementcavityD2Q9 model
格子Boltzmann方法(LBM)^[1-3]广泛应用于流体动力学计算中, 采用Boltzmann方程代替格子气自动机, 并将该模型用于流体的数值计算^[4-5], 再引入平衡态分布函数, 将碰撞算子用一个碰撞矩阵代替^[6].在此基础上, 单弛豫时间法(SRT)^[7]进一步简化了碰撞算子.Qian等^[8]提出格子BGK模型, D’humeriers^[9]提出了一种广义LBE模型(GLBE), Lallemand等^[10]对GLBE模型作了细致的理论分析, 表明其在物理原理、参数选取和数值稳定性方面都有很大优势.Luo等^[11]提出多弛豫碰撞模型在精度和数值稳定性方面都要优于单弛豫的BGK(Bhatnagar-Gross-Krok)碰撞模型, 当对高雷诺数流动进行模拟时, BGK模型会产生数据波动, 最终导致计算失稳^[12].GLBE的碰撞过程使用多个弛豫时间, 因此也称MRT(multi-relaxation-time)模型^[13-16].多参数弛豫模型是目前使用最广泛的格子Boltzmann模型, 其中Qian等^[8]提出的D2Q9模型最具有代表性, 目前已证明MRT-LBM模型不仅能提高计算稳定性, 还能提高准确性, Peng等^[17]提出浸入式边界LBM的多块模型, 对方腔流进行分块网格处理来验证正确性, 成功计算出方腔流的涡心, 结果与Ghia等^[18]的DNS计算结果值吻合.
本文通过对物理量变化剧烈的左、右上方两角区域的网格进行局部加密处理, 将流场进行分区, 各区域使用不同密度的网格划分, 同时通过区域间边界处的信息传递来实现计算的耦合, 可显著提高计算效率.
1 MRT加密算法1.1 多参数弛豫模型格子Boltzmann方法MRT-LBM方程为

(1)

式中：f为节点上的分布函数; m为矩; m^eq为矩的平衡态; M为正交转换矩阵; S为驰豫系数矩阵.
在MRT模型中, 剪切黏度和体积黏度分别为

(2)

式中: c_s²=1/3;弛豫参数S_v与S_e分别与剪切黏度和体积黏度有关.
MRT模型的主要特征是在矩空间进行碰撞, 然后在离散的速度方向进行迁移.
1.2 局部加密方法图 1为网格加密结构, 在方腔左、右上角分别设有2个加密区, 这里以左上角为例.粗网格区边界为ABC, 加密区域边界为DEF, 为了使两区域的信息得以互相传递, 在中间分别加一过渡层.加密区过渡层为ABC, 即与粗网格区边界重合; 粗网格区过渡层为GHI, 在粗网格区和加密区中, 两节点之间的点(标记为×的叉点)为插值点, 插值点上的信息可以通过相邻两点的3次样条插值^[19]得到.
图 1(Fig. 1)

图 1 网格加密结构图Fig.1 Structure sketch of local grid refinement

粗网格区c的网格尺寸为δx^c, 加密区f的网格尺寸为δx^f, 尺寸比例为n=δx^c/δx^f=δt^c/δt^f=2, 由粗细网格区域内剪切黏度和体积黏度相同可得加密区和粗网格区的弛豫参数关系为

(3)

在区域边界上的正应力与剪应力应相等, 由Chapman-Enskog分析知

(4)

根据粗细网格界面节点上对应应力分量连续, 在这2个区域中3个非平衡态e, p_xx, p_xy之间的关系为

(5)

其他3个非平衡态矩为能量平方ε, x与y方向能量通量，由文献[19-20]得

(6)

由此可知, 粗网格区与加密区之间的矩空间关系为

(7)

加密区与粗网格区碰撞后的分布函数推导见文献[20], 加密区碰撞后的分布函数为

(8)

式中:

(9)

同样, 粗网格区碰撞后的分布函数也可以由加密区碰撞后的函数得到:

(10)

式中,
2 数值计算结果与分析为了验证MRT-LBM模型进行局部加密提高计算精度的有效性, 对方腔的左、右上角2个部位的局部网格加密(DCLR), 如图 2所示, 这2个区域存在非常大的速度梯度.
图 2(Fig. 2)

图 2 二维方腔流局部网格加密(DCLR)布局Fig.2 Sketch of local grid refinement for two dimensional cavity

流动雷诺数Re=LU_w/v=1000, 其中L是方腔宽度, 这里L=N_x=N_y=128.5;U_w是上端盖的驱动速度.v为剪切黏度, 体积黏度ξ等于v, 粗网格区内的S_e^c=1.8868, S_v^c=1.8868;δt^c=1;δx^c=1;S_ε^c=1.54, S_q^c=1.9, 两个局部加密区的宽度N_xx=N_yy=39, dx^f=0.5, 由此确定的弛豫系数S_e^f= 0.641, S_v^f=0.641, 另外2个弛豫系数不变, 即S_ε^f=1.54, S_q^f=1.9.
2.1 相对速度在初始是静止状态, 经过一段充分长的时间后(超过200000个粗网格时间步长), 使得流动达到稳定.图 3为所有速度都是方腔流达到稳定状态时的结果, 将之与Ghia等^[18]的基准值和Wu等^[21]的粗网格MRT-LBM模拟结果进行对比.由图 3可知, 沿腔内垂直几何中心线和腔内水平几何中心线的瞬时速度与Ghia的基准值非常接近, 远超过Wu的模拟结果, 可知加密网格处理后的数据更加精确.
图 3(Fig. 3)

图 3 沿方腔中心线的速度Fig.3 Velocity along mid-plane

2.2 压力轮廓图如图 4所示, 左边是经过两角局部网格加密处理后的无量纲压力轮廓图, 右边是Wu所得出的无量纲压力轮廓图.通过对比, 在速度不连续的方腔内左、右上角, 粗网格处理后会存在非常多的噪声区域, 而对两角进行局部加密后的压力图则有着非常明显的降噪效果.其次, 在两角加密区和粗网格区的边界上的压力轮廓线是连续的, 这说明了本文加密边界处理的正确性和可行性.
图 4(Fig. 4)

图 4 方腔压力轮廓图Fig.4 Pressure contour of cavity (a)—DCLR; (b)—MRT-LBM.

2.3 涡量涡量分布如图 5所示.加密后的方腔无量纲涡量图清晰地表现左、右上角两个奇异角落的漩涡强度, 而粗网格所画出的无量纲涡量图在这两个角上存在较大嘈杂噪声, 两角局部网格加密可以改善涡量图整体计算精度.
图 5(Fig. 5)

图 5 方腔涡量图Fig.5 Vortex contour of cavity (a)—DCLR; (b)—MRT-LBM.

2.4 应力由式(4)可以得到τ_xx, τ_xy, 中心截面上应力如图 6所示, 应力曲线始终是连续的, 说明了应力分布图的正确性, 加密后的应力剧烈振荡的曲线明显得以缓和.
图 6(Fig. 6)

图 6 方腔沿x=13的无量纲应力Fig.6 Normalized stress of cavity along x=13

2.5 计算资源进行局部加密处理后, 网格数量从16641个加密到25767个, 增加了9126个网格.单次循环的CPU时间由4.225ms提高到了9.255ms, 牺牲了一定的计算时间.但是在内存占用方面, 程序从4164kB内存提高到了6447kB, 仅为1.55倍, 却带来了显著的计算精度提升, 这说明局部网格加密具有非常高的可行性.
3 结论1) 通过对顶盖驱动方腔的2个奇异角进行加密网格处理后, 得到的相对速度更加精确, 数值更接近Ghia的基准值, 即局部加密能更好地传递数值, 提高全场数据的精确性.
2) 从得到的压力轮廓图和涡量轮廓图知:局部网格加密能够显著改善左、右上角的轮廓分布, 噪声得到了极大改善, 从而得到更清晰更精确的压力与涡量结果.
3) 除了速度分布, 本文还对方腔内的应力分布进行了验证.对奇异角落局部网格加密能够极大减少应力振荡, 得出较为平滑的应力曲线, 在以前的研究中没有报道, 对今后的LBM研究具有一定的参考价值.
参考文献

[1]	Chen S, Doolen G. Lattice Boltzmann method for fluid flows[J].Annual Review of Fluid Mechanics, 1998, 30: 329–364.DOI:10.1146/annurev.fluid.30.1.329
[2]	Succi S. The lattice Boltzmann equation for fluid dynamics and beyond[M]. Oxford: Clarendon Press, 2001.
[3]	Yu D, Mei R, Luo L S, et al. Viscous flow computations with the method of lattice Boltzmann equation[J].Progress Aerospace Science, 2003, 39: 329–367.DOI:10.1016/S0376-0421(03)00003-4
[4]	Frish U, Hasslaeher B, Porneau Y. Lattice-gas automata for the Navier-Stokes equations[J].Physical Review Letters, 1986, 56(14): 1505–1508.DOI:10.1103/PhysRevLett.56.1505
[5]	Menamara G R, Zanetti G. Use of the Boltzmann equation to simulate lattice automata[J].Physical Review Letters, 1988, 61(20): 2332–2335.DOI:10.1103/PhysRevLett.61.2332
[6]	Higuera F J, Jimenez J. Boltzmann approach to lattice gas simulation[J].Euro-Physics Letters, 1989, 9(7): 663–668.DOI:10.1209/0295-5075/9/7/009
[7]	Chen S, Chen H D, Martinez D, et al. Lattice Boltzmann model for simulation of magneto hydrodynamics[J].Physical Review Letters, 1991, 67(27): 3776–3779.DOI:10.1103/PhysRevLett.67.3776
[8]	Qian Y H, D'Humieres D, Lallemand P. Lattice BGK models for Navier-Stokes equation[J].Euro-Physics Letters, 1992, 17(6): 479–484.DOI:10.1209/0295-5075/17/6/001
[9]	D'humieres D. Generalized lattice Boltzmann equations in AIAA rarefied gas dynamics:theory and applications[J].Progress Astronaut Aeronaut, 1992, 159: 450–458.
[10]	Lallemand P, Luo L S. Theory of the lattice Boltzmann method:dispersion, dissipation, isotropy, Galilean invariance and stability[J].Physical Review E, 2000, 61: 6546–6562.DOI:10.1103/PhysRevE.61.6546
[11]	Luo L S, Liao W, Chen X, et al. Numerics of the lattice Boltzmann method:effects of collision models on the lattice Boltzmann simulations[J].Physical Review E, 2011, 83(5): 056710.DOI:10.1103/PhysRevE.83.056710
[12]	杨帆, 施徐明, 郭雪岩, 等. 多弛豫时间格子波尔兹曼方法的分块算法[J].排灌机械工程学报, 2013(1): 56–60. ( Yang Fan, Shi Xu-ming, Guo Xue-yan, et al. Multi-block implementation of MRT-LBM[J].Journal of Drainage and Irrigation Machinery Engineering, 2013(1): 56–60.)
[13]	郭照立, 郑楚光. 格子Boltzmann方法的原理及应用[M]. 北京: 科学出版社, 2009: 46-53. ( Guo Zhao-li, Zheng Chu-guang. The principle and application of lattice Boltzmann method[M]. Beijing: Science Press, 2009: 46-53.)
[14]	Shan F, Guo X S, Tu J, et al. Multi-relaxation-time lattice Boltzmann modeling of the acoustic field generated by focused transducer[J].International Journal of Modern Physics C, 2017, 28(3): 38–52.
[15]	Humieres D, Ginzburg I, Krafczyk M, et al. Multiple relaxation-time lattice Boltzmann models in three dimensions[J].Philosophical Transaction of the Royal Society, 2002, 360: 437–451.DOI:10.1098/rsta.2001.0955
[16]	Mei R, Luo L S, Lallemand P, et al. Consistent initial conditions for lattice Boltzmann simulation[J].Computer & Fluids, 2006, 35: 855–862.
[17]	Peng Y, Shu C, Chew Y T, et al. Application of multi-block approach in the immersed boundary-lattice Boltzmann method for viscous fluid flows[J].Journal of Computational Physics, 2006, 218: 460–478.DOI:10.1016/j.jcp.2006.02.017
[18]	Ghia U, Ghia K N, Shin C T. High-resolutions for incompressible flow using the Navier-Stokes equations and a multi-grid method[J].Journal of Computational Physics, 1982, 48: 387–411.DOI:10.1016/0021-9991(82)90058-4
[19]	Farhat H, Lee J S. Fundamentals of migrating multi-block lattice Boltzmann model for immiscible mixtures in 2D geometries[J].International Journal of Multiphase Flow, 2010, 36: 769–779.DOI:10.1016/j.ijmultiphaseflow.2010.07.001
[20]	Chen S Y, Peng C, Teng Y H, et al. Improving lattice Boltzmann simulation of moving particles in a viscous flow using local grid refinement[J].Computer & Fluids, 2010, 36(10): 769–779.
[21]	Wu J S, Shao Y L. Simulation of lid-driven cavity flows by parallel lattice Boltzmann method using multi-relaxation-time scheme[J].International Journal of Numerical Method in Fluids, 2004, 46: 921–937.DOI:10.1002/(ISSN)1097-0363