引言
格子玻尔兹曼方法（Lattice Boltzmann Method, LBM）是一种介观尺度的流体数值方法：它不直接离散 Navier-Stokes 方程，而是追踪定义在规则格子上的粒子分布函数，通过”迁移—碰撞”两个简单步骤演化，宏观流场则作为分布函数的统计矩自然涌现。LBM 算法结构高度规则、天然局部、易于并行，在多孔介质流、多相流、微尺度流等领域展现出独特优势。本文从玻尔兹曼方程的 BGK 近似出发，推导 D2Q9 离散速度模型，并给出顶盖驱动方腔流的完整 Python 实现。
理论基础
从玻尔兹曼方程到 BGK 模型
介观视角下，流体由分布函数 $f(\mathbf{x}, \boldsymbol{\xi}, t)$ 描述——它给出时刻 $t$、位置 $\mathbf{x}$ 处速度为 $\boldsymbol{\xi}$ 的粒子密度。其演化遵循玻尔兹曼方程。Bhatnagar-Gross-Krook（BGK）近似把复杂的碰撞算子简化为向局部平衡态的单松弛过程：
$$\frac{\partial f}{\partial t} + \boldsymbol{\xi} \cdot \nabla f = -\frac{f - f^{\mathrm{eq}}}{\tau}$$
其中 $\tau$ 为松弛时间，$f^{\mathrm{eq}}$ 为局部 Maxwell-Boltzmann 平衡分布。对速度空间做离散，把 $\boldsymbol{\xi}$ 限制为一组有限的离散速度 $\mathbf{c}_i$，并在规则格子上沿速度方向迁移，即得到格子 BGK（LBGK）方程：
$$f_i(\mathbf{x} + \mathbf{c}_i \Delta t,, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{\Delta t}{\tau}\left[f_i - f_i^{\mathrm{eq}}\right]$$
D2Q9 模型与宏观量
二维九速（D2Q9）模型取 $\mathbf{c}_0 = (0,0)$，四个轴向速度 $(\pm 1, 0), (0, \pm 1)$ 与四个对角速度 $(\pm 1, \pm 1)$（格子单位 $c = \Delta x / \Delta t = 1$）。对应的平衡分布为：
$$f_i^{\mathrm{eq}} = w_i \rho \left[1 + 3(\mathbf{c}_i \cdot \mathbf{u}) + \frac{9}{2}(\mathbf{c}_i \cdot \mathbf{u})^2 - \frac{3}{2}|\mathbf{u}|^2\right]$$
各方向权系数如下表：
| 方向组 | 速度 $\mathbf{c}_i$ | 个数 | 权 $w_i$ |
|---|---|---|---|
| 静止 | $(0,0)$ | 1 | $4/9$ |
| 轴向 | $(\pm 1, 0), (0, \pm 1)$ | 4 | $1/9$ |
| 对角 | $(\pm 1, \pm 1)$ | 4 | $1/36$ |
宏观密度与速度由分布函数的矩给出：$\rho = \sum_i f_i$，$\rho \mathbf{u} = \sum_i \mathbf{c}_i f_i$。通过 Chapman-Enskog 多尺度展开可以证明，该离散系统在低 Mach 数极限下恢复不可压缩 Navier-Stokes 方程，且运动黏度与松弛时间满足：
$$\nu = c_s^2\left(\tau - \frac{1}{2}\right)\Delta t, \qquad c_s = \frac{c}{\sqrt{3}}$$
这一关系是 LBM 参数标定的核心：给定目标 $\Re$ 与格子分辨率，反推 $\tau$；$\tau$ 越接近 $1/2$（黏度越小、$\Re$ 越高），数值稳定性越差，需借助 MRT（多松弛时间）或熵稳定格式改善。
数值方法与实现
以下代码实现 D2Q9 方腔流：迁移用 np.roll 整体完成，侧壁与底壁采用半步长反弹格式，运动顶盖采用 Zou-He 速度边界：
import numpy as np
def lbm_cavity(nx=100, Re=1000.0, nt=12000):
"""D2Q9 顶盖驱动方腔流（BGK 单松弛模型）"""
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w = np.array([4/9] + [1/9]*4 + [1/36]*4)
U = 0.1 # 格子单位下顶盖速度（保证 Ma<0.3）
nu = U * nx / Re
tau = 3.0 * nu + 0.5
f = np.tile(w[:, None, None], (1, nx, nx)) # 初始为 rho=1 的平衡态
for _ in range(nt):
rho = np.sum(f, axis=0)
ux = np.tensordot(cx, f, axes=1) / rho
uy = np.tensordot(cy, f, axes=1) / rho
# 碰撞：BGK 松弛到局部平衡
for i in range(9):
cu = cx[i]*ux + cy[i]*uy
feq = w[i]*rho*(1 + 3*cu + 4.5*cu**2 - 1.5*(ux**2 + uy**2))
f[i] -= (f[i] - feq) / tau
# 迁移
for i in range(9):
f[i] = np.roll(np.roll(f[i], cx[i], axis=0), cy[i], axis=1)
# 左右壁面反弹
f[[1, 5, 8], 0, :] = f[[3, 7, 6], 0, :]
f[[3, 6, 7], -1, :] = f[[1, 8, 5], -1, :]
# 底壁反弹
f[[2, 5, 6], :, 0] = f[[4, 7, 8], :, 0]
# 顶盖 Zou-He 速度边界 (u=U, v=0)
rho_t = f[0, :, -1] + f[1, :, -1] + f[3, :, -1]
+ 2*(f[2, :, -1] + f[5, :, -1] + f[6, :, -1])
f[4, :, -1] = f[2, :, -1]
f[7, :, -1] = f[5, :, -1] + 0.5*(f[1, :, -1] - f[3, :, -1]) - 0.5*rho_t*U
f[8, :, -1] = f[6, :, -1] + 0.5*(f[3, :, -1] - f[1, :, -1]) + 0.5*rho_t*U
return ux, uy
ux, uy = lbm_cavity()
print(f"松弛时间 tau = {3.0*(0.1*100/1000)+0.5:.3f}")
print(f"最大速度幅值 = {np.max(np.hypot(ux, uy)):.4f} (顶盖附近应≈0.1)")由于每层迭代仅涉及相邻格点数据交换，整个算法只有碰撞（逐点局部计算）与迁移（规则数据搬移）两个核，天然适合 GPU 与大规模并行——这也是 LBM 近年在工业界快速普及的关键原因。代码中 $\tau = 0.53$ 已接近稳定边界，若继续提高 $\Re$，需改用 MRT 碰撞模型或引入湍流亚格子模型。
工程案例分析
问题：岩心孔隙尺度渗流模拟
某砂岩 CT 重构孔隙结构分辨率为 $300^3$ 体素，孔隙度 $\phi = 0.18$。用 LBM 直接模拟孔隙内单相流，计算绝对渗透率并与实验值对比。
| 参数 | 数值 | 说明 |
|---|---|---|
| 体素数 | 300³ = 2.7×10⁷ | 孔隙/固体二值标记 |
| 松弛时间 $\tau$ | 0.8 | 兼顾稳定与精度 |
| 驱动方式 | 体积力 10⁻⁵ | 避免压力边界的压缩性误差 |
| 迭代步数 | 5×10⁴ | 流量相对变化 < 10⁻⁴ 判稳 |
| 计算渗透率 | 850 mD | 实验值 780 mD，偏差 9% |
渗透率由 Darcy 定律 $k = \mu \bar{u} L / \Delta p$ 反算。结果偏差主要来自分辨率对喉道的刻画不足——喉道宽度仅 3–5 个体素时，反弹边界引起的有效孔径偏移约半个格子。采用插值反弹（interpolated bounce-back）或局部加密可显著改善。
结论与展望
本文从玻尔兹曼方程的 BGK 近似出发，建立了 D2Q9-LBGK 模型并给出可运行的方腔流代码。LBM 的精髓在于用线性迁移加局部碰撞替代非线性对流项的离散，以介观自由度的冗余换取算法结构的简洁。进阶方向包括：MRT 与级联（cascaded）碰撞模型提升高 $\Re$ 稳定性、自由能/赝势多相流模型、浸没边界法处理运动固体、以及与共轭换热模型耦合的热 LBM。对工业大算例，GPU 上的 fused kernel 与 AA-pattern 内存布局可将单卡性能推至数十亿格点更新每秒量级。
发表回复