5. 复杂几何层流流动

文档B站视频讲解：

LBM结合半反弹边界+作用力模型求解Poiseuille流动

5.1. 半反弹边界条件

大家常说格子玻尔兹曼方法（LBM）的一大优势是边界条件容易实现，其实这主要体现在半反弹边界条件的应用上——其他类型的边界条件相比传统数值方法，并没有特别明显的优势。下面我们就用通俗的方式讲讲半反弹边界条件。

半反弹边界条件是处理固体边界的一种简单又实用的方法。核心思路很直观：当粒子分布函数撞到固体边界时，沿着入射方向过来的粒子，会沿着相反方向反弹回去，反弹方向的分布函数值就用入射方向的数值来更新。

示意图如下所示：

半反弹边界示意图

具体来说，假设在某个时间步\(t\)，位置\(\boldsymbol{x}\)处的粒子分布函数\(f_\alpha (\boldsymbol{x},t)\)沿着方向\(\boldsymbol{e}_\alpha\)撞到固体边界上，那么到下一个时间步\(t+\delta_t\)，这个位置上反弹方向的分布函数就会更新为：

\[ f_{\bar{\alpha}} (\boldsymbol{x},t+\delta_t) = f_\alpha (\boldsymbol{x},t) \]

其中，\(\bar{\alpha}\)代表和\(\alpha\)完全相反的方向。

这种边界条件的好处特别明显：一是实现起来简单，计算速度快；二是能很好地保证质量守恒（不会出现粒子凭空消失或增加的情况）。最关键的是，从上面的公式能看出来，半反弹边界的计算只需要当前位置的信息，不用依赖周围网格点的数据，这让它在并行计算中效率极高。

半反弹真实壁面与计算壁面的偏差

当然，半反弹边界条件也有小缺点。比如上图里能看到，真实的壁面位置和计算中假设的壁面位置会有偏差：计算时默认壁面在两个网格点的正中间，但实际壁面可能出现在两个网格点之间的任意位置，这就会带来一点计算误差。

半反弹模拟泊肃叶流动

半反弹模拟泊肃叶流动模拟精度

上面是用半反弹边界模拟泊肃叶流动的结果，能清楚看到模拟结果和理论解贴合得非常好，而且模型依然保持二阶精度。这个算例和1.4节全反弹的算例完全一样，只是把边界条件换成了半反弹而已。从上两张图里还能发现，对于1.4节提到的泊肃叶流动，用全反弹边界时\(x\)方向的网格数是\(n+1\)个，而用半反弹边界时只需要\(n\)个。

5.2. 作用力格式

上一节我们讲了怎么用半反弹边界处理固体边界，但如果要模拟复杂多孔结构里的流动，之前用进出口压力差驱动的方式就不太方便了——因为多孔结构里很难明确划分出进出口，自然也没法准确设定进出口压力。为了解决这个问题，我们可以改用外力驱动的方式来模拟流动。

后续我们会整理LBM里所有常见的作用力模型，这里先讲最常用的一种——Guo等的作用力格式[Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E. 65(4) (2002) 046308.]。这种格式的使用很简单，只需要在LBM的演化方程里加一个源项就行：

\[f_{\alpha}^{*}(\boldsymbol{x},t)=f_{\alpha}(\boldsymbol{x},t)-\frac{1}{\tau}[f_{\alpha}(\boldsymbol{x},t)-f_{\alpha}^{\mathrm{eq}}(\boldsymbol{x},t)]+{\color[RGB]{240, 0, 0} F_{\alpha}}\]

其中，源项\(F_\alpha\)的计算公式是：

\[{\color[RGB]{240, 0, 0} F_{\alpha}=\delta _t\omega _\alpha\left( 1-\frac{1}{2\tau} \right) \left[ \frac{\boldsymbol{e}_\alpha-\boldsymbol{u}}{c_{s}^{2}}+\frac{(\boldsymbol{e}_\alpha\cdot \boldsymbol{u})\boldsymbol{e}_\alpha}{c_{s}^{4}} \right] \cdot \boldsymbol{F}}\]

这里的\(\boldsymbol{F}\)就是我们施加的外力。同时，速度的计算也要做一点小修正：

\[\rho \boldsymbol{u} = \sum_\alpha f_\alpha \boldsymbol{e_\alpha} {\color[RGB]{240, 0, 0} + \frac{\delta_t}{2} \boldsymbol{F}}\]

此时，模型对应的宏观方程就将外力项考虑进去了。

下面是用Guo作用力格式模拟泊肃叶流动的结果，能看到模拟解和理论解依然贴合得很好，模型也保持二阶精度。这个算例和5.1节的算例完全一致，边界还是半反弹边界，只是把驱动方式从“进出口压力差”换成了“外力驱动”。压力差和驱动力的对应关系很简单：\(F = \Delta P/ L\)（\(\Delta P\)是进出口压差，\(L\)是流动方向的长度）。这时泊肃叶流动的解析解可以表示为：

\[u_y(x) = \frac{F}{2\mu}x(H-x)\]

作用力驱动模拟泊肃叶流动

作用力驱动模拟泊肃叶流动模拟精度

还有个很有意思的点：用外力驱动时，流动方向的网格可以设置得非常少。比如这种泊肃叶流动，哪怕流动方向只设3个网格点，也能算出很好的结果。

5.3. 编程实现复杂几何的层流流动

半反弹边界的一大优势是它可以通过指定0表示液体，1表示固体来模拟一些复杂的多孔介质流动。我们在上一节压力驱动的泊肃叶流动的基础上开展相关的修改。主要修改如下：

1. 通过0、1数组指定多孔介质形状

主要代码如下：

shape = zeros(n, m);
% 修改成你自己的几何，这是用一个圆柱绕流表示
for i = 1:1:n
    for j = 1:1:m
        if (((i - 0.5)* dx - length/2)^2 + ((j-0.5)*dx - width/2)^2)^0.5 <= 10
            shape(i + 1,j + 1) = 1;
        end
    end
end

2. 判断需要进行边界条件处理的网格点

主要代码如下：

index = zeros(8, n, m);
for i = 1:1:8
    index(i,shape==0 & circshift(shape,[round(e(i,1)/cc),round(e(i,2)/cc)]) == 1) = 1;
end

本段代码的作用是为每个非静止的速度方向找到需要进行处理的网格点。判断方式为，如果该网格为shape值为0（表示液体），且该网格与此速度方向的反方向网格shape为1（表示固体），则需要对该网格进行处理。

3. 边界条件处理

主要代码如下：

inv_direction = [3,4,1,2,7,8,5,6];
for i = 1:1:8
    f(i,index(i,:,:) == 1) = f_old(inv_direction(i),index(i,:,:) == 1);
end

本代码的作用是对从固体中流出的速度重新赋值，赋值规则按照半反弹方式。

相关修改代码的详细讲解，请查看讲解视频。

下面展示主要结果：

圆柱绕流模拟

随机多孔介质模拟

两个平板模拟

形状JR的J模拟

还可以模拟泊肃叶流动，模型同样具有二阶精度

作用力驱动+任意几何模拟泊肃叶流动

作用力驱动+任意几何模拟泊肃叶流动模拟精度

5.4. 代码下载

5.1半反弹边界泊肃叶.zip

5.2作用力驱动泊肃叶.zip

5.3任意几何流动模拟.zip

---

如果您觉得这个项目对您有帮助，可以考虑用以下方式支持我：

• ☕ 请我喝杯咖啡

如果条件允许，欢迎捐赠支持！
每一分都是对我莫大的鼓励，让我能投入更多时间维护和更新。

• 📄 引用我的文章

如果暂时不便捐赠，适当引用我的文章也是极好的支持！
您的引用能帮助这个工作获得更多关注，同样让我感到无比欣慰。

无论哪种方式，都是对我的巨大支持！🙏感谢您让开源世界更美好！✨