格子玻尔兹曼方法(LBM)及其在微通道绕流中的应用

冯俊杰，孙 冰，姜 杰，徐 伟，石 宁

(中国石化青岛安全工程研究院化学品安全控制国家重点实验室，山东青岛 266071)

0 前言

微反应器在提高反应过程安全性、缩短反应时间、提高转化率、灵活生产等方面具有独特的优势，实现微通道流动的精确测定和控制是微反应器发挥诸多优势的保障和广泛应用的基础[1]。由于微通道内的流动具有尺度小、多尺度、相界面与边界复杂的特点，传统的计算流体力学(CFD)方法作为宏观模拟方法存在着诸多不足，而格子玻尔兹曼方法(lattice Boltzmann method，LBM)突破了传统计算方法的框架，直接从离散模型出发，通过粒子群的碰撞和迁移代替传统的连续流体模型，更接近流动的微观本质，在微流控领域具有明显的优势[2-3]。

格子玻尔兹曼方法的核心思想是将流体离散为在网格上运动的介观粒子，通过计算粒子的碰撞和迁移规律得到粒子分布函数，进而统计计算得到宏观变量如压力、速度等分布规律，创造性地实现了模拟流体运动的连续介质模型向离散模型的转变[4]。由于LBM方法基于非平衡统计物理学的Boltzmann方程，因而能成为联系微观分子尺度与宏观尺度之间的纽带[5-6]。传统的CFD方法主要基于宏观的连续介质假设，而难以计算那些不符合连续介质假设或者难以用宏观方程描述的系统，对于这些体系往往需要借助微观的分子动力学或气体动理论来进行描述[7]。对于分子动力学来说必须同时跟踪大量粒子的运动，实际求解的计算量非常大。在这种背景下，基于分子运动论和概率统计力学的LBM方法就成为一种有效方法，其具有更高的计算效率，并且容易实现并行计算[8-10]。

绕流是工程中以及自然界中常见的流体运动行为，化工反应器中的搅拌桨叶的转动、液相流过催化剂层、挡板及其他内构件等现象都属于绕流，这种广泛发生的物理现象包含着深刻的传递机理[11-12]。研究者已针对绕流的数值模拟开展了很多卓有成效的工作，有限差分、有限体积以及有限元等离散方法，大涡模拟、雷诺平均纳维-斯托克斯、离散涡等流场计算方法，被广泛应用于绕流行为的模拟研究[13-14]。

由于微通道特征尺寸较小(一般在10～1 000 μm之间)，其内部流动机理与宏观尺度存在显著的差异，宏观尺度下的许多规律不再适用于微通道中，目前对于微通道中的绕流行为理解仍不够深入，实现微反应器的流动精确控制及科学设计仍旧存在困难。

本文阐述了格子玻尔兹曼方法的原理和计算方法，并应用该方法构建数值模型进行了宏观尺度及微通道中的绕流计算模拟，结果可为相关多相流动过程提供理论指导与数据支撑。

1 LBM基本原理与方法

LBM的前身为“格子气”(Lattice Gas Automata)或“元胞自动机”模型，Boltzmann方程是统计力学中用以描述非平衡态分布函数演化规律的方程，在Boltzmann方程的求解过程中，其难点在分布函数的非线性项——碰撞项，该非线性项还与分子间的具体作用力有关。格子气自动机方法有存在随机噪音、碰撞算子的指数复杂等缺点，20世纪80年代研究者通过引入平衡态分布函数将碰撞算子线性化，降低了碰撞算子的指数复杂性，直到90年代提出的单松弛时间法和LBGK模型，进一步简化了碰撞算子，此后LBM逐渐发展成熟，近年来已成为流体力学领域研究热点之一[15]。

BGK近似通过一个简单的算子来替代碰撞项从而简化Boltzmann方程，而格子玻尔兹曼方程是Boltzmann—BGK方程的进一步离散形式，这一离散形式包括了速度离散、时间离散、空间离散。至今LBM模型中，使用最为广泛的单松弛模型或LBGK模型方程如下：

式中：τ——无量纲松弛时间；

t——时间；

fi——分布函数；

feq——局部平衡分布函数。

如果将数值粘度吸收到物理粘度中，则LBGK方程的空间和时间精度都是二阶的。对上式根据特定时间步长Δt积分可得：

LBM的粒子分布的网格模型通常用DdQm模型表达，其中前者代表网格维度，后者代表粒子可能的运动方向，平衡态分布函数可统一表示为：

式中：wi——权重因子；

c——粒子速度，m/s；

cs——格子声速，m/s。

目前，最常用的基本模型有D2Q7、D2Q9、D3Q15、D3Q19、D3Q27模型等，网格结构如图1所示。

图1 LBM方法中的主要网格模型

本文中使用的D2Q9模型的速度配置如下：

对于格子玻尔兹曼方程的初始化有两种方法——非平衡态校正方法和迭代方法，前者基本思想是基于Chapman-Enskog展开，求出分布函数的高阶项近似，进而得到分布函数的近似表达式；后者基本思路是求解初始压力的Poisson方程，同时初始化分布函数，以求得到与速度场相一致的初始分布函数。此外，边界处理是格子玻尔兹曼方法不同于传统CFD数值模拟的重要内容，根据计算格式的不同主要可分为启发式、动力学、外推以及其他复杂边界处理格式。

综上所述，完整的格子玻尔兹曼模型主要分为三个部分：格子，即离散速度模型；平衡态分布函数；分布函数的演化方程。其计算过程可以归纳为：根据研究对象特性，首先确定所有节点上的宏观物理量，进而计算各节点在各方向上的平衡态分布函数，获得初始场；求解离散后的控制方程，针对不同的边界节点特征采用不同的边界处理格式，基于格子玻尔兹曼模型的宏观量的定义法则，计算各节点的物理量分布。

2 绕流过程的LBM模拟

2.1 宏观尺度绕流模拟

本文首先对典型的宏观尺度流体绕圆柱流动行为进行模拟，采用D2Q9模型与反弹格式边界条 件(即粒子碰撞固体边界后反弹回流动域，保证边界的质量与动量守恒)。圆柱体半径为0.25 m，高度为1 m，将其置于一个x方向风速为10 m/s的风洞中。对于模拟重点关注圆柱体周围及尾流区域的格子分布进行了加密，如图2所示。

图2 圆柱体网格设置及其周围网格密度变化示意

经过模拟得到圆柱周围绕流流体速度分布(俯视图)随时间的变化如图3所示。

图3 宏观尺度绕流流场随时间变化规律

由模拟结果可知，绕流经过约2 s之后达到稳定状态，障碍物两侧周期性地脱落出旋转方向相反、并排列成有规则的双列线涡，这两列线涡分别保持自身的运动前进，进而互相干扰、吸引，形成了非线性的涡街。涡街的形成原理是由于绕流在任何时候左右都不会完全对称，某一时刻必有一方压力高，而另一方压力低，所以必定产生相对运动，与此同时，流动使得低压区也在发展，当其壮大到一定的程度后开始反向运动，至压力平衡时移动停止，导致尾涡脱落。由于流体在不断的流动，上述的过程往复进行，最终实现周期性的尾涡脱落现象。

流体力学中涡量是描写旋涡运动强度和方向的最重要的物理量之一，其数学定义为速度场的旋度，涡量的大小是流体微团绕该点旋转的平均角速度的两倍，方向与微团的瞬时转动轴线重合。本研究中绕流涡量随时间变化规律如图4所示，可以看到尾涡变化经历了对称——混沌——周期性稳定的过程，在本条件下，尾涡脱落周期约为0.1 s。

图4 绕流涡量随时间变化规律

本案例为圆柱绕流的经典体系，已有物理学家通过理论推导得出x方向上阻力系数理论值约为1，而y方向上以原点为中心呈周期性振荡[16]。通过本文模拟结果计算障碍物的受力情况，可以得到绕流经过圆柱体时x、y两个方向的阻力系数，分别如图5显示，可以发现模拟得到的x、y两个方向的阻力系数均与理论值非常接近，验证了本模拟模型的准确性。

图5 绕流阻力系数随时间变化

2.2 微通道中的绕流过程数值模拟

本节中使用2.1节建立的LBM模型，研究微通道尺度下的绕流行为，与2.1节算例保持绕流停留时间一致，圆柱体半径为0.25 mm，高度为10 mm，将其一个水平风速为0.01 m/s的风洞中。模拟得到的速度云图如图6所示。

图6 微通道绕流流场稳定后速度分布

本研究条件下微通道绕流流场与宏观尺度相比，能够迅速达到稳定状态，而且不发生周期性变化。由图可知，微通道绕流分布整体更加均匀对称，障碍物对微通道绕流的阻碍影响范围小于5倍直径范围，而且尾涡区域更接近层流的梯度分布，湍动效果不明显。计算流场的涡量分布，如图7所示。

图7 微通道绕流流场稳定后涡量分布

由涡量云图变化图可以发现，并未出现涡街脱落现象。在2.1节的案例中，由于液相流速比较高，物体两边的涡开始互相影响，也就有了所谓的卡门涡街。而本条件下液相流速比较小，物体两边产生的涡效应较弱，互相影响较小，因此难以形成周期性的涡街，最终体现为障碍物与绕流之间的作用力更加稳定。

选取圆柱中心所在水平线(x方向)计算尾流流动信息，如图8所示，其中原点位置为圆柱右侧边缘。可以发现障碍物后的速度分布呈递增规律，且增幅逐渐放缓，在障碍物后3 mm附近(约5倍直径距离)达到最大值，此后基本保持稳定。而通过涡量分布可以发现，在尾流中心线上涡量先增后减，在0.5 mm附近(约1倍直径距离)达到最大值，此后逐渐下降，同样是在障碍物后3 mm附近达到稳定值。

图8 障碍物中心线上尾流流动变化

3 结论

本文介绍的LBM方法基本的计算变量是微观上的粒子分布函数，通过求解离散玻尔兹曼方程，采用碰撞模型(如BGK模型)来模拟牛顿型流体流动，从而代替求解N-S方程，模拟得到的有限个粒子的流动和碰撞等相互作用成为宏观粘性流体流动的缩影。对绕流过程中流场结构、固体阻力、脱落周期等变化规律的分析表明，格子玻尔兹曼方法计算稳定可靠，效率高，能够用于科学地计算微反应器及微流控领域的数值模拟；同等停留时间条件下，微反应器中的圆柱绕流相对宏观尺度湍动程度明显降低，尾流未形成周期性的涡街，在实际生产过程中体现为流动更加均匀、可控，有助于实现化学反应的精确控制。