基于格子Boltzmann方法的固定方柱绕流研究

靳遵龙， 王元凯， 王永庆

(郑州大学 化工与能源学院 河南 郑州 450001)

0 引言

实际生活中很多情形都可以归结为钝体绕流问题，如河流绕过桥墩的流动以及建筑物的风载等.流体流经钝体时，由于过流界面减小，会产生一系列复杂的物理现象，特别是边界层的分离进而在钝体后产生周期脱落的涡.目前钝体绕流研究主要以圆柱绕流为主，对于一般形状柱体的研究还更多的处于实验阶段.

对于复杂流体绕流，传统的CFD方法是从宏观Navier-Stokes方程出发，利用各种数值方法进行模拟.文献[1]分别使用有限差分法和离散涡方法计算了低雷诺数和高雷诺数下的方柱绕流；文献[2]通过有限元法求解了2个圆柱左右并排时和前后排列时的绕流.格子Boltzmann方法作为一种介观模拟方法，将宏观运动和微观运动的统计平均联系起来，在流体力学问题的数值模拟方面具有广阔的前景.格子Boltzmann方法不但可以方便地模拟复杂几何形状边界的流体流动，还可以模拟系统的时间演化，并且只需要简单的算法就可以在计算机上实现并行计算，具有更高的计算效率.本文利用格子Boltzmann方法对不同雷诺数下固定单方柱绕流流场进行分析，并对并列双方柱不同分布间距的流场进行了模拟，验证了格子Boltzmann方法边界处理和数值模拟的正确性和便捷性.

1 格子Boltzmann方法

图1 D2Q9模型离散速度Fig.1 Discrete velocities of D2Q9 model

格子Boltzmann方法中的基本变量是格点的分布函数fi，fi为沿i方向的粒子分布函数.本文采用D2Q9模型，9个离散速度分布如图1所示.其演化方程(不含外力项)为

式中：Δt为时间步长；τ为无量纲松弛时间，τ=τ0/Δt；离散速度ci表示粒子的运动方向，可以表示为

c0=(0,0)，c1,3,c2,4=(±c,0),(0,±c)，c5,6,c7,8=(±c,c),(c,±c)，

由平衡态分布函数通过Chapman Enskog展开，可以导出模型的宏观密度和宏观速度分别为

格子Boltzmann方法在物理空间上将系统粒子运动分为迁移和碰撞两个相对独立的过程，使其具备很好的并行特性以及较强的复杂边界处理能力.

2 格子Boltzmann方法的初始条件和边界格式处理

流体的运动总是在一定的初始条件和边界条件下进行的，因此必须设置相应的初始条件和边界条件.格子Boltzmann方法中的基本变量是分布函数.对于稳态或者准稳态流动，初始条件对最终计算结果影响不大，可以直接将初始分布函数设为其平衡态分布函数.实际问题中，边界条件往往基于宏观物理量，如何根据宏观量合理地确定分布函数是格子Boltzmann方法中的重要问题.常用的边界主要分为启发格式边界、动力学格式边界、外推格式边界以及复杂格式边界；根据实际问题又可分为速度边界和压力边界等.

2.1 非平衡态外推格式边界

2.2 非平衡态反弹格式边界

实际应用中，边界处的速度分量或者压力分布往往是已知的，如管道流动等.文献[4]在1997年提出的非平衡态反弹格式则可以很好地应用于上述两种边界条件. 对于速度边界条件，有

通过非平衡态反弹格式可以直接计算垂直于速度边界上的3个未知平衡态分布函数，相似的方法也可以用于压力边界和其他边界条件.

3 固定方柱绕流的数值模拟与分析

3.1 单方柱绕流模拟

图2 单方柱绕流模型示意图Fig.2 Schematic diagram of channel flow past a single square cylinder

图3为不同雷诺数下流场涡量云图.可以看出，当Re较小时，流场是定常的，方柱后面出现一对上下对称的尾涡.随着Re的增加，当Re=38时，流场逐渐转变为非定常，尾涡也逐渐脱落，尾流区中形成周期性摆动和交错的漩涡，即出现卡门涡街.随着Re的进一步增加，当Re=238时，尾涡有紊乱的趋势.进一步的模拟显示，当Re大于238时，尾流区由层流状态过渡到紊流状态.

图3 不同雷诺数下流场涡量云图Fig.3 Vorticity contours at different Reynolds numbers

图4显示了Re=200时阻力系数和升力系数随时间的变化曲线，当前平均阻力系数和最大升力系数分别为1.51和0.31，与文献[5]中实验研究的数据相比具有一致性，可见利用格子Boltzmann方法模拟是可行的.

图4 Re=200时阻力系数和升力系数随时间的变化曲线Fig.4 Time history of drag coefficient and lift coefficient at Re=200

Re=35及Re=37时的流场流线分析如图5所示，可以看出，上下对称的尾涡逐渐扩大，上涡有分离的趋势.因此，定常流失稳的临界Re在37左右，这与文献[6]中的结论是相符合的.

图5 不同雷诺数下流场流线图Fig.5 Streamline patterns at different Reynolds numbers

与理论研究及文献中圆柱绕流的结果对比后发现，方柱绕流整体上符合一般的钝体绕流趋势，并且由于方柱的特殊性，临界Re又与圆柱绕流不同.由以上分析可知，方柱绕流中卡门涡街现象在Re为40～238时较为明显.本文将采用Re=200进行下一步的模拟研究.

3.2 双方柱绕流模拟

图6 并列双方柱绕流模型示意图Fig.6 Schematic diagram of channel flow past parallel square cylinders

实际生活中更多的是流体绕过多个柱体流动的情形，柱体的个数、形状、排列方式等均会影响流场的结构.本文研究双方柱并列模型下不同分布间距对流场的影响，控制的变量为柱间距M与方柱特征尺寸H的比值(M/H)，并列双方柱绕流模型示意图如图6所示.

采用相同的边界条件，X、Y方向格子数分别为360和120，在Re=200的条件下分别对并列方柱M/H=1、M/H=2、M/H=3、M/H=4四种条件进行模拟，其涡量云图如图7所示.

图8显示了并列方柱在M/H=2和M/H=3 时的流场流线图.从图7和图8中可以看出，当分布间距较小(M/H=1)时，两柱形成的漩涡彼此影响较小，还保留着部分单柱绕流的流场特征；随着分布间距的增加，流场逐渐变得复杂，漩涡分布呈现不规则状态，当M/H=2时，流体绕流方柱的涡流彼此影响最为明显；分布间距继续增加，当M/H大于3时，两个涡的相互影响越来越小，方柱下游形成两个反向同步脱落的涡街，各自涡形逐渐接近单方柱卡门涡街状态.

图7 不同分布间距下流场涡量云图Fig.7 Vorticity contours at different distribution distances

图8 不同分布间距下流场流线图Fig.8 Streamline patterns at different distribution distances

4 小结

通过对单方柱和并列双方柱的模拟研究，验证了格子Boltzmann方法边界处理和数值模拟的正确性和便捷性.通过对不同Re下单方柱流场的模拟分析，探究了方柱产生卡门涡街的临界雷诺数(Re=37)，并得到了产生明显卡门涡街的雷诺数范围(Re=40～238).此外，雷诺数Re=200下并列双方柱不同分布间距(M/H=1、M/H=2、M/H=3、M/H=4)的流场模拟结果表明，当柱间距为2倍方柱边长时，流体绕流方柱的涡流彼此影响最为明显.