格子Boltzmann方法三种边界格式的对比分析

格子Boltzmann方法三种边界格式的对比分析

刘连国;杨帆;王宏光

【摘要】采用格子Boltzmann方法(Lattice Boltzmann Method- LBM)对二维

顶盖驱动方腔流动进行数值模拟.在计算中分别使用半步长反弹、非平衡反弹、以

及非平衡外推三种边界处理格式,并得到了不同格式对应的流线分布,流函数最小值、涡心坐标、几何中心线速度分布等.通过将所得结果与基准解进行比较,就三种边界

格式的计算效率,计算精度、以及计算稳定性等方面进行了讨论和分析,为LBM计

算中边界格式的选择提供了有益的参考.

【期刊名称】《机械研究与应用》

【年(卷),期】2012(000)001

【总页数】5页(P18-22)

【关键词】格子Boltzmann方法;边界处理格式;半步长反弹格式;非平衡反弹格式;非平衡外推格式

【作者】刘连国;杨帆;王宏光

【作者单位】上海理工大学能源与动力工程学院,上海200093;上海理工大学能源

与动力工程学院,上海200093;上海理工大学能源与动力工程学院,上海200093【正文语种】中文

【中图分类】O357.1

1 引言

格子Boltzmann方法(LBM)是近年来迅速发展的一种新型数值计算方法。边界条

件的处理是LBM实施中一项非常关键的内容。实际计算表明:选取不同的边界条件会对数值计算的精度、稳定性以及效率产生很大影响。作为LBM的一个基本问题，边界条件的处理一直是流体力学一个重要的研究方面。根据边界条件的类型，可将之分为两类:压力边界和速度边界［1］，其中的速度边界又可细分为:平直边界和

曲面边界。笔者从经典的流体力学问题二维顶盖方腔流模拟入手，对三种平直边界格式进行对比和分析，为LBM计算中边界格式的选择提供了有益的参考。

2 二维九点格子Boltzmann模型

目前最常用的格子Boltzmann模型为LBGK模型，通过引入“单一弛豫时间”来简化Boltzmann方程中碰撞项的计算［2］。

九点格子LBGK模型的演化方程为:

式中:(x，t)是在t时刻、x处的平衡态分布函数;τ为单一弛豫时间因子;eα为网格

点各方向上的粒子速度。

其中:c=δx/δt。每个格点处的密度ρ和宏观速度u满足如下方程:

压力由状态方程求得:

其中:声速c s=c/。该模型中，运动粘度v与弛豫时间τ的关系为:

v=(τ-1/2)粒子平衡态分布函数为［3］:

3 三种边界格式简介

3.1 标准反弹格式

对于静止固体边界，常用的处理方法为对边界上的粒子作弹回处理，称为标准反弹格式。

图1 标准反弹格式

如图1所示，自流体节点(i-1，2)入射(迁移)到边界节点(i，1)的分布函数f8，不发生碰撞即沿原路弹回，由此可获得节点(i，1)上的f7。类似地，边界节点(i，1)的

其他未知分布函数也可由相应流体节点的分布函数经过反弹获得，即:

式中:f4，6，8(i，1)由 f4(i，2)，f6(i+1，2)和 f8(i-1，2)迁移得到。

标准反弹格式操作简单，能严格保持系统的质量和动量守恒，却只有一阶精度［4］，降低了格子Boltzmann方法的精度。

3.2 半步长反弹格式

为改进标准反弹格式的精度，1993年，Ziegler提出了半步长反弹格式［5］。该格式引入了计算边界的概念。该格式与标准反弹格式完全相同，也是在紧邻边界的流体格点上执行碰撞，但是物理边界不是放置在格点上，而是位于计算网格的中心，如图2所示，半步长反弹格式具有清晰的物理图像:在(i，2)处碰撞后，速度指向壁面的粒子经过δt/2后到达壁面，与壁面碰撞并反转其速度，再经δt/2到达格点(i，2)，因此有:

图2 半步长反弹格式

图3 非平衡反弹格式

3.3 非平衡反弹格式

1997年，Zou与He对二维和三维的LBGK模型进行研究，提出一种新的边界处理格式，称为非平衡反弹格式(nonequilibrium bounce-back scheme)［6］。

如图3所示，2-0-1位于边界上，6-4-8和7-3-5分别位于非流体区域(如固体壁面)内和流体区域内。当每一个迁移步完成后，在边界节点O处，分布函数f0、f1、f2、f4、f6 和 f8 已知，而 f3、f5、f7 和密度ρ未知。以速度边界为例，流体宏

观速度的水平分量ux和垂直分量uy已知，可得如下方程组:

以上有4个未知数，但只有3个独立方程。为使方程封闭，Zou和He假设与边

界垂直的方向上，反弹格式对分布函数的非平衡态部分仍然成立，即:

联立方程式(5)～(8)，可求得边界节点O处的所有未知量。非平衡反弹格式可取得近似二阶精度的模拟结果。

3.4 非平衡外推格式

Guo等人于2002年提出了一种新的边界处理格式，即非平衡外推格式(nonequilibrium extrapolation scheme)［7］，其基本思想是，将边界节点上

的分布函数分解为平衡态和非平衡态两部分，其中平衡态部分由边界条件的定义近似获得，而非平衡态部分则用非平衡外推确定。

图4 非平衡外推格式

如图4所示，假设COA位于边界上，EBD位于流场中。在每次碰撞之前，需要知道边界节点O上的分布函数fα(O，t)，并将其分为平衡态和非平衡态两部分，即:

对于平衡态部分(O，t)，可用边界节点上的宏观变量求得，如果节点O存在未知

宏观变量，则由B点的相应值代替。以速度边界条件为例，已知O点速度u(O，t)，密度ρ(O，t)未知，则O点的平衡态分布函数由下式近似获得:

对于非平衡态部分(O，t)，由于流体节点 B处的分布函数fα(O，t)、宏观速度