基于LBM方法的方柱绕流数值计算

基于Boltzmann方法的方柱绕流数值模拟

张华

（武汉理工大学 1049721101590）

摘要：基于格子Boltzmann方法( Lattice Boltzmann Method，LBM) ，本文对在较高Re数下的方柱绕流问题进行了数值模拟。给出了正方形柱体在一定攻角下的流线图。在较低Re数下数值结果较好，高Re数下，数值计算呈现不稳定性。

关键字：格子Boltzmann方法；方柱绕流；数值模拟； D2Q9模型

Boltzmann Method Based On Numerical Simulation Of

Fluid Flows Around Square Cylinder

Hua Zhang

（WHUT 1049721101590）

Abstract: Based on the lattice boltzmann method (LBM), this paper simulated the problem of fluid flows around squares in higher Re number. And it is showed stream line chart of the square in a certain angle of attack. As a result, the numerical results in lower Re number are more better then that in higher Re number.

Key words: lattice boltzmann method； Flow around a square cylinder； numerical simulation； D2Q9 model

1 引言

模拟流体运动的CFD数值计算方法主要有两种:(l)从宏观角度出发,基于连续介质假设,采用数值计算方法,求解Euler方程或N-S方程;(2)从微观角度出发,采用分子动力学的方法,对流动进行数值模拟[1]。

传统的计算流体力学中有限差分、有限元、有限体积等数值计算方法都属于第一种方法,这种方法直接对非线性的微分方程例如熟知的Euler方程和N-S方程进行离散,得到代数方程组或微分方程系统,然后在用标准的数值方程求解。虽然近些年CFD取得了很大的进展，但由于流体运动的复杂性和计算机资源的限

制，这类方法比较直观但也存在许多不足。

第二种数值方法不同于第一种，它是一种全新的数值模拟方法。它采用分子动力学的方法，从微观的角度对流动进行数值模拟。这种方法就是近来发展的格子波尔兹曼方法( Lattice BoltzmannMethod ，LBM) ，它是一种基于粒子分布函数演化的数值方法。在格子波尔兹曼方法中，由于粒子速度与空间离散一致，粒子在空间网格点上传输，对流仅仅是简单的赋值运算，没有空间插值误差。粒子分布函数通过粒子在网格点的碰撞而改变，仅由当地网格点上宏观值决定，不涉及其它网格点的值，因而具有优秀的局部性，便于并行计算。目前，LBM 已经用于多个领域［2］，如湍流、多向流、粒子流、磁流体力学、多孔介质渗透流动以及一些具有复杂边界的流动[3][4]。

2 格子Boltzmann 方法的基本原理

LBM(Lattice Boltzmann,格子Boltzmann 法) 作为一种数值方法，它不同于传统的有限差分法、有限元法及有限体积法，在LBM 中，基本的计算变量不是密度、速度这些宏观的物理量，而是细观上的粒子分布函数。通过粒子分布函数的演化分布，从细观层次上研究流动的宏观现象。格子Boltzmann 方法是一种在空间和时间上离散，介于连续与离散之间的数值方法。从格子 Boltzmann 方法出发，可以在一定的条件下推导出 Navier-Stokes 方程。

2.1 控制方程

一旦求得粒子分布函数,则宏观速度和压强就可以由其前的两个动量自动求出。LBM 中粒子分布函数),(t x f a ,表示在时刻t ,位置x ,速度为αe 的粒子分布,其演化可以分解为以下两步(BGK 模型)[5][6] :

碰撞: τ

δαααα),(),(),(),(~)(t x f t x f t x f t t x f eq -+=+ （1） 迁移：

),(~),(t t x f t t t e x f δδδααα+=++ （2）

式中，τ为松弛时间，),()(t x f eq α是平衡函数。宏观的密度、动量由下式求得

∑∑∑∑

====)()(eq eq f e f e u f f α

αααααα

αρρ （3）

粘性系数为 t c ∆-=)5.0(3

2

τν （4） 其中，t x c ∆∆=/为粒子速度，声速3/c C S =。这里的τ可取2

1Re 32+∆=∞t c D U τ 2.2 二维9速度模型（D2Q9）[6][7]

LB 方法建立模型的核心问题是根据网格形式确定与之对应的平衡分布函数表达式。而只有当平衡分布函数已知的情况下才能进一步演化LB 方法的动力方程。一般情况下，LB 方法各种模型所剖分的网格具有物理对称性，其中包括了平衡分布函数中权重组合的对称和各个参数的选择。因此，不同的网格剖分形式有着不同的平衡分布函数。

图1是正方形网格二维9 速度(D2Q9，D 指维数，Q 指粒子运动方向的总数)模型。整个流场剖分为正方形网格，每个节点与周围 8 个节点相邻，加上零

速度，粒子共有9个运动方向，所有9个方向上的速度矢量构成一个集合αe 。

图1 二维 9 速度(D2Q9)模型

根据速率的大小，将αe 分成三类：

1）();0,0,0==ααe

2）;4,3,2,1,21sin ,21cos =⎪⎭⎫ ⎝⎛--=απαπααe

3）;8,7,6,5,425sin ,425cos 2=⎪⎪⎭⎫ ⎝⎛⎪⎭⎫ ⎝⎛+-⎪⎭⎫ ⎝⎛+-=αππαππααe

在不可压缩条件下，利用D2Q9模型完全可以模拟 N-S 方程。得到不可压缩条件下的D2Q9 模型的平衡分布函数：

8,7,6,5,23)(29)(313614,3,2,1,23)(29)(31910

,2319422)(22)(2)(=⎥⎦

⎤⎢⎣⎡-⋅+⋅+==⎥⎦

⎤⎢⎣⎡-⋅+⋅+==⎥⎦⎤⎢⎣⎡-=αραραραααααααu u e u e f u u e u e f u f eq eq eq （5） 将式（5）写成统一形式：

⎥⎦⎤⎢⎣

⎡-⋅+⋅+=22)(23)(29)(31u u e u e f eq ααααρω （6） 其中

⎪⎩⎪⎨⎧====.

8,7,6,53614,3,2,19104αααωα，；，

；， 格子 Boltzmann 方法是一种与传统数值方法不同的流体计算和建模方法。它最大的优点是运用简单的线性运算模拟各种复杂的非线性宏观现象。状态方程就是压力和密度成一个线性关系，不需要另外的压力修正，压力表示如下

ρ2s c p = （7）

2.3 边界处理方法

边界条件的处理方法在格子Bolztmann 方法中起重要作用,对格子Bolztmann 模型的精度和稳定性都有很大的影响。这里我们介绍两种文中用到的边界条件处理方法：周期格式和反弹格式[6]。