格子玻尔兹曼方法(LBM)模拟液滴在重力下穿孔（相场模型）C++代码

最近在研究多相流相关的内容，尝试使用格子玻尔兹曼方法（LBM）来模拟液滴在重力作用下穿孔的过程，这里采用相场模型来处理。不得不说，这个过程就像一场充满挑战与惊喜的编程冒险，今天就来和大家分享一下其中的关键代码和实现思路。

基本原理速览

格子玻尔兹曼方法基于微观粒子的运动和碰撞来描述宏观流体的行为。在相场模型里，我们通过引入一个相场变量来区分不同的相，比如液滴和周围的介质。而重力的引入则为整个模拟增添了现实世界的物理特性，让液滴有了“掉落”和“穿孔”的动力。

C++ 代码解析

初始化部分

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

const int Lx = 200; // 模拟区域在x方向的长度
const int Ly = 200; // 模拟区域在y方向的长度
const double tau = 1.0; // 松弛时间
const double rho0 = 1.0; // 初始密度
const double g = 0.01; // 重力加速度

std::vector<std::vector<std::vector<double>>> f(Lx, std::vector<std::vector<double>>(Ly, std::vector<double>(9, 0.0)));
std::vector<std::vector<double>> rho(Lx, std::vector<double>(Ly, rho0));
std::vector<std::vector<std::vector<int>>> e = {
    {0, 0}, {1, 0}, {0, 1}, {-1, 0}, {0, -1}, {1, 1}, {-1, 1}, {-1, -1}, {1, -1}
};

在这段代码里，我们首先包含了必要的头文件，为后续的输入输出、向量操作和数学运算做准备。然后定义了模拟区域的大小（Lx 和 Ly），这就像是搭建了一个虚拟的“舞台”，液滴的各种行为都将在这里上演。松弛时间 tau 是格子玻尔兹曼方法里的一个关键参数，它决定了分布函数向平衡态松弛的速度。初始密度 rho0 和重力加速度 g 也都设定好了，这俩参数分别给了液滴初始的“物质含量”和“掉落的动力”。

格子玻尔兹曼方法(LBM)模拟液滴在重力下穿孔（相场模型）C++代码

我们用三维向量 f 来存储每个格点上不同方向的分布函数，二维向量 rho 存储每个格点的密度。还有 e 向量，它定义了不同方向的离散速度，这在后续的粒子运动计算中至关重要。

平衡分布函数计算

double equilibrium(int i, double rho_val, double ux, double uy) {
    double cu = 3.0 * (e[i][0] * ux + e[i][1] * uy);
    double u2 = 3.0 * (ux * ux + uy * uy);
    if (i == 0) {
        return 4.0 / 9.0 * rho_val * (1.0 - 0.5 * u2);
    } else if (i >= 1 && i <= 4) {
        return 1.0 / 9.0 * rho_val * (1.0 + cu + 0.5 * cu * cu - 0.5 * u2);
    } else {
        return 1.0 / 36.0 * rho_val * (1.0 + cu + 0.5 * cu * cu - 0.5 * u2);
    }
}

平衡分布函数是格子玻尔兹曼方法的核心之一。这个函数根据给定的格点方向 i、密度 rho_val 以及速度分量 ux 和 uy 来计算平衡态下的分布函数值。不同的方向 i 对应不同的计算公式，大致思路是基于宏观的密度和速度信息，推算微观层面上粒子在各个方向的分布情况。

碰撞步骤

void collision() {
    for (int x = 0; x < Lx; ++x) {
        for (int y = 0; y < Ly; ++y) {
            double ux = 0.0, uy = 0.0;
            for (int i = 0; i < 9; ++i) {
                ux += e[i][0] * f[x][y][i];
                uy += e[i][1] * f[x][y][i];
            }
            ux /= rho[x][y];
            uy /= rho[x][y];

            for (int i = 0; i < 9; ++i) {
                double feq = equilibrium(i, rho[x][y], ux, uy);
                f[x][y][i] = (1.0 - 1.0 / tau) * f[x][y][i] + 1.0 / tau * feq;
            }
        }
    }
}

碰撞步骤模拟了粒子在格点上的相互作用。首先，我们在每个格点计算宏观速度 ux 和 uy，方法是对各个方向的分布函数加权求和再除以密度。然后，根据当前格点的密度和速度计算每个方向的平衡分布函数 feq，最后通过松弛过程，让当前的分布函数 f 向平衡态靠近，这就完成了一次碰撞模拟。

流步骤

void streaming() {
    std::vector<std::vector<std::vector<double>>> f_new(Lx, std::vector<std::vector<double>>(Ly, std::vector<double>(9, 0.0)));
    for (int x = 0; x < Lx; ++x) {
        for (int y = 0; y < Ly; ++y) {
            for (int i = 0; i < 9; ++i) {
                int x_new = (x + e[i][0] + Lx) % Lx;
                int y_new = (y + e[i][1] + Ly) % Ly;
                f_new[x_new][y_new][i] = f[x][y][i];
            }
        }
    }
    f = f_new;
}

流步骤描述了粒子从一个格点移动到相邻格点的过程。我们先创建一个新的分布函数数组 fnew，用来存储移动后的分布函数。对于每个格点和每个方向 i，根据离散速度 e[i] 计算粒子移动后的新位置 xnew 和 y_new，并把原格点的分布函数值赋给新格点。最后更新 f，让模拟继续向前推进。

重力处理

void apply_gravity() {
    for (int x = 0; x < Lx; ++x) {
        for (int y = 0; y < Ly; ++y) {
            double ux = 0.0, uy = 0.0;
            for (int i = 0; i < 9; ++i) {
                ux += e[i][0] * f[x][y][i];
                uy += e[i][1] * f[x][y][i];
            }
            ux /= rho[x][y];
            uy /= rho[x][y];

            uy += g;

            for (int i = 0; i < 9; ++i) {
                double cu = 3.0 * (e[i][0] * ux + e[i][1] * uy);
                double u2 = 3.0 * (ux * ux + uy * uy);
                double feq;
                if (i == 0) {
                    feq = 4.0 / 9.0 * rho[x][y] * (1.0 - 0.5 * u2);
                } else if (i >= 1 && i <= 4) {
                    feq = 1.0 / 9.0 * rho[x][y] * (1.0 + cu + 0.5 * cu * cu - 0.5 * u2);
                } else {
                    feq = 1.0 / 36.0 * rho[x][y] * (1.0 + cu + 0.5 * cu * cu - 0.5 * u2);
                }
                f[x][y][i] += 1.0 / tau * (feq - f[x][y][i]);
            }
        }
    }
}

在这个函数里，我们对每个格点重新计算速度，然后给垂直方向速度分量 uy 加上重力加速度 g。接着重新计算平衡分布函数，并基于新的平衡态对分布函数进行松弛，这样就把重力的影响引入到了模拟中，促使液滴开始“掉落”。

总结与展望

通过上述代码和分析，我们完成了使用格子玻尔兹曼方法（LBM）结合相场模型模拟液滴在重力下穿孔的初步实现。当然，这只是一个基础版本，还有很多可以优化和拓展的地方，比如边界条件的精细处理、相场模型参数的进一步调优，让液滴的行为更加符合实际物理现象。希望这篇博文能给同样在探索这个领域的小伙伴一些启发，一起在代码的世界里解锁更多有趣的物理模拟。