【LBM教程 1】格子 Boltzmann 方法概述

关注微信同名公众号，获取更多信息

前言

格子 Boltzmann 方法（Lattice Boltzmann Method，简称 LBM）作为数值流体力学（CFD）领域中一种新兴而高效的模拟方法，近年来受到广泛关注与应用。从传统的有限差分、有限体积等方法到基于分子运动统计理论的 LBM，其独特的理论基础和数值实现方式为复杂流动问题的研究提供了全新的视角和工具。本文作为本系列教程的第一篇文章，旨在为读者构建对 LBM 方法的初步认识，从理论渊源、基本概念、核心思想以及部分代码示例等多个角度进行详细介绍，为后续各篇深入探讨奠定坚实基础。

在传统的 CFD 数值方法中，研究者通常需要求解 Navier–Stokes 方程，这是一组高度非线性的偏微分方程。尽管这些方法在许多工程应用中取得了成功，但在处理复杂边界条件、多相流、湍流以及微尺度流动时，常常面临求解效率低、稳定性差等问题。相比之下，LBM 采用了介于微观粒子运动和宏观连续介质描述之间的中间模型，通过离散化的分子分布函数来刻画流体的统计行为，其数值格式简单、易于并行计算，能够更好地适应大规模复杂几何和多物理耦合问题的模拟需求。

LBM 的理论基础与基本思想

格子 Boltzmann 方法的起源可以追溯到对 Boltzmann 方程的离散化。Boltzmann 方程描述了气体分子在碰撞与迁移过程中的概率分布演化，其数学表达式为:
$$\frac{\partial f}{\partial t}+\mathbf{v}\cdot \nabla f+\mathbf{F}\cdot {\nabla }_{\mathbf{v}}f=\Omega (f)$$
其中 $f=f(\mathbf{x},\mathbf{u},t)$ 表示在位置 $\mathbf{x}$、速度 $\mathbf{u}$ 和时刻 $t$ 的粒子分布函数，$\mathbf{F}$ 为外力作用，而 则表示碰撞项。直接求解 Boltzmann 方程非常复杂，因此 LBM 通过对速度空间和物理空间进行离散化，将连续的分布函数转化为离散的“格子”分布函数 $f_i(\mathbf{x},t)$，其中下标 $i$ 对应于预先确定的离散速度方向。

最常见的离散模型之一是 D2Q9 模型，其在二维平面内构造 9 个离散速度方向；类似的，在三维问题中常用 D3Q19 或 D3Q27 模型。这种离散化方法使得碰撞和迁移过程都可以在格子点上以局部操作的形式进行，从而大大简化了数值实现过程，并且天然适合于并行计算。

在 LBM 中，碰撞过程通常采用 BGK（Bhatnagar-Gross-Krook）碰撞模型进行近似。BGK 模型的核心思想是，粒子分布函数在碰撞过程中以一个单一的松弛时间 $\tau$ 向平衡分布函数 $f_i^{eq}$ 收敛，其离散形式可以写为：
$$f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)=f_i(\mathbf{x},t)-\frac{1}{\tau }\left[f_i(\mathbf{x},t)-f_i^{eq}(\mathbf{x},t)\right]$$

其中 ${\mathbf{c}}_i$ 为离散速度，$\Delta t$ 是时间步长。平衡分布函数 $f_i^{eq}$ 通常通过对麦克斯韦分布函数做适当展开得到，其形式能够保证质量和动量等守恒定律在碰撞过程中的正确实现。

从上式中可以看出，LBM 的计算主要分为两个步骤：碰撞和迁移。在碰撞步骤中，每个格点上的分布函数根据局部信息与平衡状态进行松弛，完成微观状态的更新；而在迁移步骤中，各个格点的分布函数沿着预定的离散速度方向传递到相邻的格点，从而实现信息的空间传播。这种局部化操作不仅使得算法简单高效，还使得对边界条件的处理变得更为灵活。

数值离散化与 LBM 的优越性

传统方法求解 Navier–Stokes 方程时，常常需要借助大量的数值求解技巧来平衡稳定性和精度，而 LBM 则利用了离散速度模型将流体运动分解成简单的局部更新过程。例如，在 D2Q9 模型中，每个格点存储 9 个分布函数，且碰撞与迁移操作只涉及局部的数据交换，从而极大地降低了算法复杂度和编程难度。与此同时，由于每个格点的计算过程都是独立的，LBM 自然具备良好的并行计算性能，非常适合在 GPU 或多核 CPU 环境下实现大规模数值模拟。

LBM 的另一个显著优势在于其处理复杂边界条件的灵活性。在传统 CFD 方法中，复杂边界条件常常需要通过复杂的网格生成与插值技术来解决，而在 LBM 中，边界处理通常通过简单的“反弹法（Bounce-Back）”或“Zou-He 边界条件”等方法实现。例如，在反弹法中，当分布函数遇到固体边界时，可以直接将其沿入射方向反弹，从而实现无滑移边界条件的模拟。这种方法不仅直观简单，而且在许多实际应用中都表现出了良好的效果。

平衡分布函数的构造

平衡分布函数 $f_i^{eq}$ 在 LBM 中起着至关重要的作用。为了保证宏观物理量（如密度 $\rho$ 和速度 $\mathbf{u}$）的正确恢复，通常要求满足以下关系：
$$\rho =\sum _i{f}_i,\rho \mathbf{u}=\sum _i{\mathbf{c}}_i{f}_i$$

在实际推导中，平衡分布函数可以通过对连续麦克斯韦分布函数进行低阶矩展开得到。在 D2Q9 模型下，平衡分布函数的经典表达式为:
$$f_i^{eq}=w_i\rho \left[1+\frac{{\mathbf{c}}_i\cdot \mathbf{u}}{c_s^2}+\frac{({\mathbf{c}}_i\cdot \mathbf{u}{)}^2}{2c_s^4}-\frac{\mathbf{u}\cdot \mathbf{u}}{2c_s^2}\right]$$

其中 $w_i$ 是各方向的权重系数，$c_s$ 为声速。这个公式不仅能够保证局部守恒定律，同时在低速流动情况下能够较好地近似连续流体的物理行为。
在碰撞阶段，BGK 模型以一个松弛时间 $\tau$ 控制分布函数向平衡状态收敛的速率，公式中 $\tau$ 的取值直接影响数值稳定性和物理粘性系数 $\nu$ 的关系，通常有:
$$\nu =c_s^2\left(\tau -\frac{1}{2}\right)\Delta t$$
这一关系不仅揭示了碰撞过程与宏观粘性的内在联系，也为后续参数选择和稳定性分析提供了理论依据。

LBM 的优势与实际应用前景

从理论到实现，LBM 展现出许多传统 CFD 方法难以比拟的优点。首先，LBM 基于离散化的分布函数描述流体运动，不再直接求解偏微分方程，极大地降低了算法复杂度；其次，碰撞和迁移过程均为局部运算，这不仅使得程序代码简洁，而且天然适合并行化处理，可以利用 GPU 或多核 CPU 进行高效计算；此外，LBM 在处理复杂边界条件、多相流、湍流、热传导等问题时，提供了灵活且有效的数值方法，已在环境科学、生物医学、能源工程等领域中得到广泛应用。

实际工程应用中，LBM 常用于模拟微流体、孔隙介质流动、流固耦合以及流体中的化学反应等问题。其在处理复杂几何边界、非均质介质以及动态变化系统中的表现尤为出色。例如，在多相流模拟中，LBM 能够自然地捕捉界面演化，而在湍流模拟中，通过适当的湍流模型扩展，LBM 同样能够反映出复杂流场的细节。正因如此，越来越多的科研人员和工程师开始将 LBM 应用于实际工程问题的数值仿真，并取得了令人瞩目的成果。

C++ 示例代码初探

为了帮助读者更直观地理解 LBM 的实现过程，下面给出一个简单的 C++ 代码片段，展示如何初始化一个二维格子并进行一次碰撞和迁移操作。代码中以 D2Q9 模型为例，简化了部分边界处理和数据结构设计，目的是说明基本思路。

#include <iostream>
#include <vector>
#include <cmath>

// 定义网格尺寸
const int NX = 100;
const int NY = 100;

// 离散速度个数（D2Q9模型）
const int Q = 9;

// 离散速度方向及对应权重
int cx[Q] = {0, 1, 0, -1, 0, 1, -1, -1, 1};
int cy[Q] = {0, 0, 1, 0, -1, 1, 1, -1, -1};
double w[Q] = {4.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0,
               1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0};

double tau = 0.6; // 松弛时间

// 初始化分布函数数组
std::vector<std::vector<std::vector<double>>> f(NX, std::vector<std::vector<double>>(NY, std::vector<double>(Q, 0.0)));

// 计算平衡分布函数
double feq(int k, double rho, double ux, double uy) 
{
    double cu = cx[k]*ux + cy[k]*uy;
    double usqr = ux*ux + uy*uy;
    double cs2 = 1.0/3.0;
    return w[k] * rho * (1.0 + 3.0*cu + 4.5*cu*cu - 1.5*usqr);
}

// 主函数
int main() 
{
    // 初始化均匀流场：密度为1，初始速度为0
    double rho0 = 1.0;
    double ux0 = 0.0, uy0 = 0.0;
    for (int i = 0; i < NX; i++) 
    {
        for (int j = 0; j < NY; j++) 
        {
            for (int k = 0; k < Q; k++) 
            {
                f[i][j][k] = feq(k, rho0, ux0, uy0);
            }
        }
    }
    
    // 一次碰撞和迁移操作示例
    // 创建临时数组保存碰撞后的分布函数
    auto f_temp = f;
    for (int i = 0; i < NX; i++) 
    {
        for (int j = 0; j < NY; j++) 
        {
            // 计算宏观变量
            double rho = 0.0, ux = 0.0, uy = 0.0;
            for (int k = 0; k < Q; k++) 
            {
                rho += f[i][j][k];
                ux += f[i][j][k] * cx[k];
                uy += f[i][j][k] * cy[k];
            }
            ux /= rho; uy /= rho;
            
            // 碰撞步：松弛向平衡分布函数
            for (int k = 0; k < Q; k++) 
            {
                double feq_val = feq(k, rho, ux, uy);
                f_temp[i][j][k] = f[i][j][k] - (f[i][j][k] - feq_val) / tau;
            }
        }
    }
    
    // Streaming 步：将碰撞后的分布函数传递到相邻格点
    std::vector<std::vector<std::vector<double>>> f_next = f;
    for (int i = 0; i < NX; i++) 
    {
        for (int j = 0; j < NY; j++) 
        {
            for (int k = 0; k < Q; k++) 
            {
                int in = (i + cx[k] + NX) % NX;
                int jn = (j + cy[k] + NY) % NY;
                f_next[in][jn][k] = f_temp[i][j][k];
            }
        }
    }
    
    // 输出某一点的密度值作为简单验证
    int test_i = NX/2, test_j = NY/2;
    double test_rho = 0.0;
    for (int k = 0; k < Q; k++) 
    {
        test_rho += f_next[test_i][test_j][k];
    }
    std::cout << "中心点密度: " << test_rho << std::endl;
    
    return 0;
}

上述代码展示了如何在二维网格上初始化分布函数、计算平衡分布函数、以及利用 BGK 模型进行碰撞和迁移的基本流程。虽然其中省略了复杂的边界处理和并行化设计，但这段示例代码足以帮助读者理解 LBM 模拟的基本步骤和算法框架。