不可压缩绕流的并行格子Boltzmann模拟

张建影, 闫广武

(1. 长春工业大学 基础科学学院, 长春 130012； 2. 吉林大学 数学学院, 长春 130012)

0 引 言

格子Boltzmann方法(LBM)是一种模拟不可压缩黏性流的工具, 已得到广泛关注[1-2]. 与传统计算流体动力学(CFD)的方法不同, LBM从细观尺度上的Maxwell方程出发, 通过特殊的离散方式给出粒子束的演化方程, 进一步得到宏观物理量. 因此, 最初的对流扩散方程, 如Navier-Stokes方程, 被一个格子Boltzmann方程取代. 由于标准格子Boltzmann方法具有如下优点： 1) 由状态方程给出压力, 不需要求解压力方程[3]； 2) 公式简单易于实现编程； 3) 算法容易实现并行[4]； 4) 复杂边界处理较方便[5]. 因此, LBM在流体力学和非线性偏微分方程等领域应用广泛[6]. 但当流动区域较大或需要处理流动细节的计算问题时, 计算网格数量的要求已超出普通计算机的能力. 在这种情况下, 需要使用并行模式的格子Boltzmann 算法[7-8], 利用LBM与其他算法的结合实现并行计算[9-12].

本文提出一个简单的并行模式, 通过在流动区域内划分若干子块以及与周围子块之间的重叠边界, 每个子块的计算是并行的[13], 利用MPI(message passing interface)并行平台实现全场的并行计算, 进而达到增加计算网格数量的目的.

1 格子Boltzmann 模型

在标准的二维网格空间中, 定义fα(x,t)为位置x、 在时刻t沿方向α的粒子分布函数. 流场中的密度和动量定义为

(5)

平衡态分布函数为

(6)

(7)

其中τ为松弛因子.

引入Knudsen数ε, 其定义为平均自由程l与特征长度L之比, 即

(8)

假设ε与格子尺度Δx同量级, 并令Δt=ε, 则Δx=c·Δt=c·ε, 利用Chapman-Enskog展开及时间多尺度展开, 可得不同时间尺度下的系列偏微分方程[14]：

将式(9)+ε×式(10)并求和得

在不可压缩假设和低速度假设下, 可得不可压缩流动的Navier-Stokes方程:

其中

(15)

2 管道内绕流的并行格子Boltzmann算法

图1 矩形计算区域子块划分Fig.1 Sub-block partition of rectangular computing area

将矩形区域用D2Q9网格铺满, 记水平方向为x轴m个格点, 竖直方向为y轴n个格点. 为了采用MPI并行技术, 将计算区域分解成如图1所示的若干子块, 并要求一个子块与相邻子块之间有重合的两列网格, 第一个子块和最后一个子块只有一个相邻的子块, 子块的数量根据全场y方向网格数决定. 在计算中, 这些子块同时计算, 除第一个和最后一个子块外, 每块数据都需要调用5部分, 分别为左右两列以及该块的内部区域[13].

计算图1中第一个子块B1时, 需要调用的数据为： 列M1N1,M2N2,M5N5,M6N6和块M3N3N4M4. 时刻(t+Δt)完成数据M2N2N5M5的计算, 列M1N1和M5N5作为该子块的边界使用. 考虑到子块之间的对接, 采用双边界M5N5,M6N6, 即在计算第一个子块时M6N6是边界, 但在计算第二个子块B2时M6N6为内部, 此时边界为M5N5. 第二个子块B2的计算需要调用列M5N5,M6N6,M9N9,M10N10和块M7N7N8M8. 列M5N5,M10N10为边界.

2.1 方柱绕流

流场为600×100网格, 运算分为3个子块, 每块为200×100. 用区域[92,108]×[42,58]表示边长16个网格的方柱, 中心位于(100,50)处, 上下边界采用固体边界. 图2为Re=100, 时间步分别为t=12 000,18 000,24 000,30 000时的涡线； 图3为Re=300, 时间步分别为t=12 000,18 000,24 000,30 000时的涡线.

图2 Re=100时方柱绕流涡量等值线Fig.2 Vorticity isolines of flow around a square cylinder at Re=100

由图2可见, 当Re=100时, 方柱后面开始出现Karmann涡街, 流动变为不定常, 受其影响, 方柱所受的升力出现周期振荡, 方柱的阻力也出现振荡. 由图3可见, 当Re=300时, 方柱后面的涡旋区域变大, 流动很快变成不定常. 这些结果与已有的数值结果[15-16]相符.

图3 Re=300时方柱绕流涡量等值线Fig.3 Vorticity isolines of flow around a square cylinder at Re=300

阻力系数与升力系数分别定义为

(16)

其中:FD为阻力;FL为升力;ρ∞为来流密度;U∞为来流速度;d为方柱特征宽度.

图4为Re=100时阻力系数与升力系数随时间步的变化曲线. 由图4可见, 当Re=100时, 阻力系数约为1.31, 而升力系数约为0, 与文献[15-16]中的阻力系数和升力系数相同, 当t>20 000时出现振荡. 图5为Re=300时阻力系数与升力系数随时间步的变化曲线. 此时方柱后面Karmann涡街的区域变大, 由图5可见， 当t>13 000时, 升力与阻力的振荡振幅变大, 与文献[15-16]的结果吻合.

图4 Re=100时方柱绕流阻力系数和 升力系数与时间步的关系曲线Fig.4 Relationship curves between drag coefficient and lift coefficient of flow around a square cylinder and time step at Re=100

图5 Re=300时方柱绕流阻力系数和 升力系数与时间步的关系曲线 Fig.5 Relationship curves between drag coefficient and lift coefficient of flow around a square cylinder and time step at Re=300

2.2 圆柱绕流

流场为600×100网格, 运算分为3个子块, 每块为200×100. 圆柱中心位于(100,50), 半径为12.5, 上下边界采用固体边界. 图6为Re=300, 时间步分别为t= 20 000,30 000,40 000,50 000时的涡线. 由图6可见, Karmann涡街中涡旋较密集. 图7为Re=300时圆柱绕流的阻力系数与升力系数随时间步的变化曲线. 由图7可见, 升力和阻力均出现振荡, 与文献[17]中的LBM结果相符.

2.3 NACA0012机翼绕流

流场为600×300网格, 运算分为3个子块, 每块为200×300. 机翼前缘位于(150,150), 翼弦长为75个网格, 攻角12°, 上下边界采用Neumann边界. 图8为Re=1 000时圆柱绕流的阻力系数与升力系数随时间步的变化曲线. 由图8可见, 升力和阻力均出现振荡, 与文献[18-19]的结果相符. 图9为Re=1 000, 时间步分别为t=6 000,8 000,10 000,40 000时的涡线. 由图9可见, 当t>6 000时机翼后面出现涡旋, 当t>10 000时出现较稳定的Karmann涡街.

图6 Re=300时圆柱绕流涡量等值线Fig.6 Vorticity isolines of flow around a circular cylinder at Re=300

图7 Re=300时圆柱绕流阻力系数和 升力系数与时间步的关系曲线Fig.7 Relationship curves between drag coefficient and lift coefficient of flow around a circular cylinder and time step at Re=300

图8 Re=1 000时圆柱绕流阻力系数和 升力系数与时间步的关系曲线 Fig.8 Relationship curves between drag coefficient and lift coefficient of flow around a circular cylinder and time step at Re=1 000

图9 Re=1 000, 攻角α=12°时涡量等值线Fig.9 Vorticity isolines at Re=1 000 and attack angle α=12°

综上, 本文应用并行模式的格子Boltzmann算法模拟了方柱、 圆柱以及NACA0012机翼绕流问题, 并给出了不同雷诺数下不同时间的涡线图, 由这些涡线图可观察到被绕物体后面形成的Karmann涡街, 由升力阻力系数随时间的变化曲线可观察到升力阻力均出现振荡, 这些结果都与经典的数值结果相符， 显示了本文方法的有效性.

[1] CHEN Shiyi, Doolen G D. Lattice Boltzmann Method for Fluid Flows [J]. Annual Review of Fluid Mechanics, 1998, 30: 329-364.

[2] Aidun C K, Clausen J R. Lattice-Boltzmann Method for Complex Flows [J]. Annual Review of Fluid Mechanics, 2010, 42: 439-472.

[3] HE Xiaoyi, CHEN Shiyi, Doolen G D. A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit [J]. J Comput Phys, 1998, 146(1): 282-300.

[4] Schmieschek S, Shamardin L, Frijters S, et al. LB3D: A Parallel Implementation of the Lattice-Boltzmann Method for Simulation of Interacting Amphiphilic Fluids [J]. Computer Physics Communications, 2017, 217: 149-161.

[5] Rosis A, de. Ground-Induced Lift Enhancement in a Tandem of Symmetric Flapping Wings: Lattice Boltzmann-Immersed Boundary Simulations [J]. Comput Struct, 2015, 153: 230-238.

[6] ZHANG Jianying, YAN Guangwu. Lattice Boltzmann Model for Complex Ginzburg-Landau Equation in Curvilinear Coordinates [J]. Computers and Mathematics with Applications, 2015, 70(12): 2904-2919.

[7] Clausen J R, Reasor D A, Jr, Aidun C K. Parallel Performance of a Lattice-Boltzmann/Finite Element Cellular Blood Flow Solver on the IBM Blue Gene/P Architecture [J]. Comput Phys Commun, 2010, 181(6): 1013-1020.

[8] GAO Jinfang, XING Huilin, Rudolph V, et al. Parallel Lattice Boltzmann Computing and Applications in Core Sample Feature Evaluation [J]. Transport in Porous Media, 2015, 107(1): 65-77.

[9] Campos J O, Oliveira R S, Santos R W, dos， et al. Lattice Boltzmann Method for Parallel Simulations of Cardiac Electrophysiology Using GPUs [J]. Journal of Computational and Applied Mathematics, 2016, 295: 70-82.

[10] Zudrop J, Masilamani K, Roller S, et al. A Robust Lattice Boltzmann Method for Parallel Simulations of Multicomponent Flows in Complex Geometries [J]. Computers and Fluids, 2017, 153: 20-33.

[11] Randles A, Kaxiras E. Parallel in Time Approximation of the Lattice Boltzmann Method for Laminar Flows [J]. Journal of Computational Physics, 2014, 270: 577-586.

[12] SU Yan, Ng T, Davidson J H. A Parallel Non-dimensional Lattice Boltzmann Method for Fluid Flow and Heat Transfer with Solid-Liquid Phase Change [J]. International Journal of Heat and Mass Transfer, 2017, 106: 503-517.

[13] 闫广武. 二维直角弯道黏性流动的格子气体仿真 [D]. 长春: 吉林大学, 1988. (YAN Guangwu. Lattice Gas Simulation of 2-D Viscous Flow in a Right-Angled Bend [D]. Changchun: Jilin University, 1988.)

[14] 闫广武. 二维剪切流的格子Boltzmann模拟 [J]. 水动力学研究与进展, 2003, 18(5): 521-525. (YAN Guangwu. Lattice Boltzmann Simulation of the Two-Dimensional Shear Flows [J]. Journal of Hydrodynamics, 2003, 18(5): 521-525.)

[15] Breuer M, Bernsdorf J, Zeiser T, et al. Accurate Computations of the Laminar Flow Past a Square Cylinder Based on Two Different Methods: Lattice-Boltzmann and Finite-Volume [J]. International Journal of Heat and Fluid Flow, 2000, 21: 186-196.

[16] Jafari S, Salmanzadeh M, Rahnama M, et al. Investigation of Particle Dispersion and Deposition in a Channel with a Square Cylinder Obstruction Using the Lattice Boltzmann Method [J]. Journal of Aerosol Science, 2010, 41(2): 198-206.

[17] Shu C, Liu N, Chew Y T. A Novel Immersed Boundary Velocity Correction-Lattice Boltzmann Method and Its Application to Simulate Flow Past a Circular Cylinder [J]. Journal of Computational Physics, 2007, 226: 1607-1622.

[18] 邱全辉, 王健平. 不可压缩机翼绕流的有限谱法计算 [J]. 计算力学学报, 2005, 22(5): 518-523. (QIU Quanhui, WANG Jianping. Finite Spectral Computation of Incompressible Flows around Airfoil [J]. Chinese Journal of Computational Mechanics, 2005, 22(5): 518-523.)

[19] Ashraf M A, Young J, Lai J C S. Reynolds Number, Thickness and Camber Effects on Flapping Airfoil Propulsion [J]. Journal of Fluids and Structures, 2017, 27(2): 145-160.