C++实现LBM模拟Couette流

使用C++实现格子玻尔兹曼方法(LBM)模拟Couette流动的完整代码，包含周期边界条件和曲线边界条件处理。

1. 基本LBM框架

#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <array>using namespace std;// 常数定义
const int NX = 100;          // x方向格子数
const int NY = 50;           // y方向格子数
const int Q = 9;             // D2Q9模型
const double tau = 0.8;      // 松弛时间
const double rho0 = 1.0;     // 初始密度
const double u0 = 0.1;       // 顶板速度// D2Q9模型的离散速度
const array<int, Q> cx = {0, 1, 0, -1,  0, 1, -1, -1,  1};
const array<int, Q> cy = {0, 0, 1,  0, -1, 1,  1, -1, -1};// 权重系数
const array<double, Q> w = {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
};// 分布函数
vector<vector<array<double, Q>>> f(NX, vector<array<double, Q>>(NY));
vector<vector<array<double, Q>>> f_eq(NX, vector<array<double, Q>>(NY));
vector<vector<array<double, Q>>> f_new(NX, vector<array<double, Q>>(NY));// 宏观量
vector<vector<double>> rho(NX, vector<double>(NY, rho0));
vector<vector<double>> ux(NX, vector<double>(NY, 0.0));
vector<vector<double>> uy(NX, vector<double>(NY, 0.0));// 边界条件标志 (0=流体, 1=上壁面, 2=下壁面, 3=曲线边界)
vector<vector<int>> boundary(NX, vector<int>(NY, 0));// 初始化函数
void initialize() {// 设置边界条件for (int i = 0; i < NX; ++i) {boundary[i][0] = 2;       // 下壁面 (静止)boundary[i][NY-1] = 1;    // 上壁面 (移动)}// 初始化分布函数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] = w[k] * rho[i][j];}}}
}// 计算平衡分布函数
void calculate_eq() {for (int i = 0; i < NX; ++i) {for (int j = 0; j < NY; ++j) {double u2 = ux[i][j] * ux[i][j] + uy[i][j] * uy[i][j];for (int k = 0; k < Q; ++k) {double cu = cx[k] * ux[i][j] + cy[k] * uy[i][j];f_eq[i][j][k] = w[k] * rho[i][j] * (1.0 + 3.0 * cu + 4.5 * cu * cu - 1.5 * u2);}}}
}// 碰撞步骤
void collision() {for (int i = 0; i < NX; ++i) {for (int j = 0; j < NY; ++j) {if (boundary[i][j] == 0) { // 仅处理流体节点for (int k = 0; k < Q; ++k) {f_new[i][j][k] = f[i][j][k] - (f[i][j][k] - f_eq[i][j][k]) / tau;}}}}
}

2. 边界条件处理

// 周期边界条件处理
void apply_periodic_boundary() {// 左边界 (i=0)for (int j = 0; j < NY; ++j) {if (boundary[0][j] == 0) { // 流体节点for (int k = 0; k < Q; ++k) {if (cx[k] == -1) { // 向左移动的粒子f_new[NX-1][j][k] = f_new[0][j][k];}}}}// 右边界 (i=NX-1)for (int j = 0; j < NY; ++j) {if (boundary[NX-1][j] == 0) { // 流体节点for (int k = 0; k < Q; ++k) {if (cx[k] == 1) { // 向右移动的粒子f_new[0][j][k] = f_new[NX-1][j][k];}}}}
}// 曲线边界条件处理 (反弹格式)
void apply_curved_boundary() {for (int i = 0; i < NX; ++i) {for (int j = 0; j < NY; ++j) {if (boundary[i][j] == 3) { // 曲线边界节点// 找到最近的流体节点int i_fluid = i, j_fluid = j;if (i > 0 && boundary[i-1][j] == 0) i_fluid = i-1;else if (i < NX-1 && boundary[i+1][j] == 0) i_fluid = i+1;if (j > 0 && boundary[i][j-1] == 0) j_fluid = j-1;else if (j < NY-1 && boundary[i][j+1] == 0) j_fluid = j+1;// 反弹格式for (int k = 0; k < Q; ++k) {int k_opp = (k + Q/2) % Q; // 相反方向f_new[i][j][k] = f_new[i_fluid][j_fluid][k_opp];}}}}
}// 移动壁面边界条件 (Couette流上壁面)
void apply_moving_wall() {for (int i = 0; i < NX; ++i) {// 上壁面 (j=NY-1)if (boundary[i][NY-1] == 1) {ux[i][NY-1] = u0; // 设置壁面速度uy[i][NY-1] = 0.0;// 计算壁面处的平衡分布double u2 = ux[i][NY-1] * ux[i][NY-1] + uy[i][NY-1] * uy[i][NY-1];for (int k = 0; k < Q; ++k) {double cu = cx[k] * ux[i][NY-1] + cy[k] * uy[i][NY-1];f_new[i][NY-1][k] = w[k] * rho[i][NY-1] * (1.0 + 3.0 * cu + 4.5 * cu * cu - 1.5 * u2);}}// 下壁面 (j=0)if (boundary[i][0] == 2) {ux[i][0] = 0.0; // 静止壁面uy[i][0] = 0.0;// 计算壁面处的平衡分布double u2 = ux[i][0] * ux[i][0] + uy[i][0] * uy[i][0];for (int k = 0; k < Q; ++k) {double cu = cx[k] * ux[i][0] + cy[k] * uy[i][0];f_new[i][0][k] = w[k] * rho[i][0] * (1.0 + 3.0 * cu + 4.5 * cu * cu - 1.5 * u2);}}}
}

3. 流动步骤和宏观量计算

// 流动步骤
void streaming() {for (int i = 0; i < NX; ++i) {for (int j = 0; j < NY; ++j) {if (boundary[i][j] == 0) { // 仅处理流体节点for (int k = 0; k < Q; ++k) {int i_next = (i + cx[k] + NX) % NX; // 周期边界处理int j_next = j + cy[k];if (j_next >= 0 && j_next < NY) {if (boundary[i_next][j_next] == 0) { // 流向流体节点f[i_next][j_next][k] = f_new[i][j][k];}}}}}}
}// 计算宏观量
void calculate_macro() {for (int i = 0; i < NX; ++i) {for (int j = 0; j < NY; ++j) {if (boundary[i][j] == 0) { // 仅处理流体节点double rho_temp = 0.0;double ux_temp = 0.0;double uy_temp = 0.0;for (int k = 0; k < Q; ++k) {rho_temp += f[i][j][k];ux_temp += f[i][j][k] * cx[k];uy_temp += f[i][j][k] * cy[k];}rho[i][j] = rho_temp;ux[i][j] = ux_temp / rho_temp;uy[i][j] = uy_temp / rho_temp;}}}
}

4. 主程序和输出

// 添加曲线边界 (示例: 圆形障碍物)
void add_curved_boundary() {int center_x = NX/2;int center_y = NY/2;int radius = NY/4;for (int i = 0; i < NX; ++i) {for (int j = 0; j < NY; ++j) {int dx = i - center_x;int dy = j - center_y;if (dx*dx + dy*dy <= radius*radius) {boundary[i][j] = 3; // 曲线边界}}}
}// 输出结果到文件
void output(int step) {ofstream file("velocity_" + to_string(step) + ".dat");file << "VARIABLES = \"X\", \"Y\", \"U\", \"V\", \"Rho\"\n";file << "ZONE I=" << NX << ", J=" << NY << ", F=POINT\n";for (int j = 0; j < NY; ++j) {for (int i = 0; i < NX; ++i) {file << i << " " << j << " " << ux[i][j] << " " << uy[i][j] << " " << rho[i][j] << "\n";}}file.close();
}int main() {// 初始化initialize();add_curved_boundary(); // 添加曲线边界// 主循环int max_step = 10000;int output_interval = 1000;for (int step = 0; step < max_step; ++step) {calculate_eq();        // 计算平衡分布collision();           // 碰撞步骤apply_periodic_boundary(); // 应用周期边界apply_curved_boundary();  // 应用曲线边界apply_moving_wall();   // 应用移动壁面边界streaming();           // 流动步骤calculate_macro();     // 计算宏观量// 输出结果if (step % output_interval == 0) {output(step);cout << "Step: " << step << endl;}}return 0;
}

5. 代码说明和参数调整

5.1 关键参数说明

1. NX, NY: 计算域大小
2. tau: 松弛时间，影响流体粘性 (ν = (τ-0.5)/3)
3. u0: 上壁面移动速度，决定剪切率
4. w: D2Q9模型的权重系数
5. cx, cy: D2Q9模型的离散速度方向

5.2 边界条件实现

1. 周期边界: 在x方向实现周期性流动
2. 曲线边界: 使用反弹格式处理复杂边界
3. 移动壁面: 上壁面以速度u0移动，下壁面静止

5.3 扩展功能

1. 要模拟不同几何形状，修改add_curved_boundary()函数
2. 要改变流动条件，调整壁面速度或初始条件
3. 要输出更多物理量，扩展output()函数

5.4 编译和运行

使用C++编译器编译并运行：

g++ lbm_couette.cpp -o lbm_couette -O3
./lbm_couette

输出文件可以使用Tecplot、Paraview等软件可视化。

参考代码 用C++模拟LBM的couette流源代码 www.youwenfan.com/contentcsl/56505.html

这个实现提供了Couette流动的完整LBM模拟框架，包含周期边界和曲线边界处理，可以作为更复杂LBM模拟的基础。