基于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中粒子分布函数),(txfa,表示在时刻t,位置x,速度为e的粒子分布,其演化可以分解为以下两步(BGK模型)[5][6] :

碰撞:

),(),(),(),(~)(txftxftxfttxfeq （1） 迁移： ),(~),(ttxftttexf （2） 式中，为松弛时间，),()(txfeq是平衡函数。宏观的密度、动量由下式求得 )()(eqeqfefeuff







（3）

粘性系数为 tc)5.0(3

2 （4）

其中，txc/为粒子速度，声速3/cCS。这里的可取21Re32tcDU 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,425cos2e 在不可压缩条件下，利用D2Q9模型完全可以模拟 N-S 方程。得到不可压缩条件下的D2Q9 模型的平衡分布函数：

8,7,6,5,23)(29)(313614,3,2,1,23)(29)(31910,2319422)(22)(2)(uueuefuueuefuf

eqeqeq



（5）

将式（5）写成统一形式： 

22)(23)(29)(31uueuefeq

 （6）

其中







.8,7,6,53614,3,2,191094，；，；，

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

2scp （7）

2.3 边界处理方法 边界条件的处理方法在格子Bolztmann方法中起重要作用,对格子Bolztmann模型的精度和稳定性都有很大的影响。这里我们介绍两种文中用到的边界条件处理方法：周期格式和反弹格式[6]。 2.3.1 周期格式 周期格式是指，当流体粒子从一个边界离开流场时，在下一时间步就从流场的另一侧进入流场。周期性边界条件主要用于空间上周期变化的流场或者无穷大的流场。周期边界条件可以表示如下：

),(),(0txfttxfNii （8） 其中0x为空间上最初始节点的位置，Nx为空间上最末节点的位置。if表示碰撞后的分布函数。即节点Nx碰撞后的分布函数就是作为节点0x下一时间步碰撞前的分布函数。

2.3.2 反弹格式 反弹格式主要用于无滑移边界条件的处理，即物面条件。当一个流体节点上的粒子沿粒子运动方向流动到下一步后到达边界节点，则下一步该粒子沿原方向的相反方向反弹回原流体节点。在这里，固壁节点不参与碰撞。反弹格式的实现如图2。

图2 反弹格式示意图 A点是流体内部的节点。A点所处的下一层节点为固壁节点，流过后，格点A的五个方向的分布函数87431,,,,fffff可由临近的格点E,D,F,C,B确定。

根据壁面反弹思想，有 867542,,ffffff . （9）

如果允许在边界格点上也进行碰撞，所得到的格式称为修正反弹格式。在该格式中，碰撞前边界上的未知分布函数设置为反方向的分布函数，即 ),(),(txftxfii （10）

处理无滑移边界的另一种改进反弹格式是半步长反弹格式。该格式中边界不是放置在格点上，而是置于格线中点上。在靠近壁面的第一层格点上仍然执行标准的反弹格式。本文中采用标准反弹格式。

3 LBM数值模拟 3.1 LBM数值模拟的简要步骤 粒子的运动只需用迁移和碰撞两种简单运算就可描述，其中它的碰撞过程是一种比较单一的时间松弛过程，粒子以单一的松弛时间逐渐趋于平衡态。即它的计算过程是一个不断循环演化的过程，与此同时，边界条件的处理在时间和空间上有很大的灵活性，压力可以直接用代数进行求解，整个算法过程是简单的显式迭代。

运用格子Boltzmann方法对流畅进行数值模拟的步骤可以归纳如下： （1） 首先进对其行空间离散，并确定相应的格子模型； （2） 给每个节点的分布函数赋以初值，从而计算它的速度和压力； （3） 根据确定的格子模型和宏观（速度及压力）物理量，计算该时刻节点的平衡分布函数；

（4） 在每一个时刻每个节点的所有方向上，根据 LB 方程对分布函数演进;再采取适当的边界处理方法求出边界点上未知方向的分布函数；

（5） 依据所得到下一时刻每个节点上的分布函数，重复以后的过程，直到满足所需所有的计算条件。

3.2 方柱绕流数值模拟 3.2.1 计算区域模型 如下图3所示30D×8D的矩形区域，将方柱放置流场中间左端为入口边界，右端为出口边界，方柱左侧边界距流场入口为8D，方柱距离下端3.5D ，D为方柱直径，网格划分300 × 80。