用格子玻尔兹曼方法 - 浸没边界法模拟圆柱绕流（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 模拟方法有了更深入的认识也欢迎一起探讨更多相关的技术细节。