用格子玻尔兹曼方法 - 浸没边界法模拟圆柱绕流（LBM - IBM in C++）
格子玻尔兹曼方法-浸没边界法模拟圆柱绕流 LBM- IBM （C++）
在计算流体力学（CFD）的领域里，格子玻尔兹曼方法（Lattice Boltzmann Method, LBM）和浸没边界法（Immersed Boundary Method, IBM）的结合为模拟复杂边界条件下的流动现象提供了一种强大的手段。今天咱就来讲讲如何用 C++ 实现基于 LBM - IBM 的圆柱绕流模拟。
格子玻尔兹曼方法基础
LBM 是一种介观尺度的数值方法，它基于对流体分子运动的简化描述。在 LBM 中，流体被看作是由大量在离散格子上运动的粒子组成。这些粒子在格子间按照特定的规则迁移和碰撞。
常见的二维格子模型是 D2Q9 模型，这里“D2”表示二维空间，“Q9”表示每个格子点有 9 个离散的速度方向。下面是 D2Q9 模型速度方向的定义代码：
// 定义 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};这里cx和cy数组分别存储了每个速度方向在 x 和 y 方向上的分量。
在 LBM 中，核心的方程是格子玻尔兹曼方程（LBE），它描述了粒子分布函数在时间和空间上的演化。简化后的离散形式如下：
\[ fi(x + ci \Delta t, t + \Delta t) - fi(x, t) = -\frac{1}{\tau} (fi(x, t) - f_i^{eq}(x, t)) \]
格子玻尔兹曼方法-浸没边界法模拟圆柱绕流 LBM- IBM （C++）
其中 \( fi \) 是在位置 \( x \)、时间 \( t \) 沿速度方向 \( ci \) 的粒子分布函数，\( f_i^{eq} \) 是平衡态分布函数，\( \tau \) 是松弛时间。
平衡态分布函数 \( f_i^{eq} \) 可以通过下面的公式计算：
\[ fi^{eq}(\rho, u) = wi \rho \left( 1 + \frac{ci \cdot u}{cs^2} + \frac{(ci \cdot u)^2}{2 cs^4} - \frac{u^2}{2 c_s^2} \right) \]
这里 \( wi \) 是权重系数，\( cs \) 是格子声速，一般在 D2Q9 模型中 \( cs = \frac{1}{\sqrt{3}} \)。权重系数 \( wi \) 的定义如下：
// D2Q9 模型的权重系数 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 radius = 10.0; // 圆柱半径 double center_x = domain_size_x / 2.0; double center_y = domain_size_y / 2.0; // 圆柱中心位置在每一步计算中，我们需要判断哪些格子点位于圆柱内部，哪些在外部。这可以通过计算格子点到圆柱中心的距离来实现：
for (int i = 0; i < nx; ++i) { for (int j = 0; j < ny; ++j) { double dx = i - center_x; double dy = j - center_y; double distance = std::sqrt(dx * dx + dy * dy); if (distance < radius) { // 该格子点在圆柱内部 // 后续可在此处施加浸没边界条件相关操作 } } }圆柱绕流模拟的整体流程
- 初始化：设置计算域大小、时间步长、流体初始条件等，初始化粒子分布函数 \( f_i \)。
- 碰撞步骤：根据 LBE 计算平衡态分布函数 \( fi^{eq} \)，并更新粒子分布函数 \( fi \)，考虑碰撞项。
- 迁移步骤：将粒子分布函数按照速度方向迁移到相邻的格子点。
- 施加浸没边界条件：在圆柱边界处施加虚拟力，调整粒子分布函数。
- 计算宏观物理量：如密度 \( \rho \) 和速度 \( u \)。
- 重复步骤 2 - 5直到达到设定的时间步数。
下面是一个简化的主循环代码框架：
// 主循环 for (int time_step = 0; time_step < total_time_steps; ++time_step) { // 碰撞步骤 for (int i = 0; i < nx; ++i) { for (int j = 0; j < ny; ++j) { // 计算平衡态分布函数 // 更新粒子分布函数 } } // 迁移步骤 for (int i = 0; i < nx; ++i) { for (int j = 0; j < ny; ++j) { // 按照速度方向迁移粒子分布函数 } } // 施加浸没边界条件 for (int i = 0; i < nx; ++i) { for (int j = 0; j < ny; ++j) { // 判断格子点是否在圆柱边界 // 施加虚拟力 } } // 计算宏观物理量 for (int i = 0; i < nx; ++i) { for (int j = 0; j < ny; ++j) { // 计算密度和速度 } } }通过上述步骤和代码实现，我们就可以用 LBM - IBM 方法在 C++ 中模拟圆柱绕流现象了。当然，实际的代码实现还需要考虑更多细节，如边界条件处理、数值稳定性等。但这个基本框架为理解和进一步开发相关模拟程序提供了一个良好的起点。希望大家对这个有趣的 CFD 模拟方法有了更深入的认识，也欢迎一起探讨更多相关的技术细节。