简介这是一份基于格子玻尔兹曼方法LBM的三维渗流模拟MATLAB源码主要面向流体力学、多孔介质渗流方向的工程师与科研人员适用于地下水流动、石油开采、土壤污染扩散等典型渗流场景。该方法通过离散玻尔兹曼方程模拟微观粒子运动进而恢复宏观流动行为对复杂边界和非均匀介质有天然优势三维算法相比二维模拟能呈现更全面的流场信息更贴近工程实际。压缩包为zip格式共1个文件为m类型脚本整体仅3KB属于轻量级核心代码方便直接阅读算法流程与参数配置。目前已有601人学习下载。借助该源码使用者可自定义流体参数、网格尺寸与边界条件执行初始化分布函数、碰撞、迁移、边界处理等标准步骤最终导出流速、压力分布的可视化流场图并计算无量纲渗流系数为渗流特性分析提供重要数据支撑。代码结构简洁既适合作为三维LBM渗流教学演示也可在此基础上做科研二次开发。1. 三维LBM渗流孔隙尺度的流动模拟入门到底难在哪“LBM三维”和“LBM渗流”连在一起指向的是一类很具体的需求用格子玻尔兹曼方法Lattice Boltzmann MethodLBM在孔隙尺度上把多孔介质中的渗流过程算出来。跟达西尺度的黑匣子模型不同三维LBM能还原每个孔隙里的速度场和压力场天然适配岩芯、页岩、土壤、多孔电极这类非均质结构。对刚入门的人它是一套可以完全用NumPy搭起来的数值框架对从二维转向三维的工程师最大的成本不在算法而在模型选型、边界处理和稳定性调参这三件事。这篇文章按我实际跑通的经验把D3Q19模型、多孔介质生成、反弹边界到渗透率输出这条链路拆开讲代码可以直接抄。2. 三维LBM选型D3Q19、松弛时间与稳定性边界2.1 三维速度模型怎么选D3Q15、D3Q19、D3Q27的取舍三维LBM比二维多出来的第一个选择就是离散速度模型。D3Q15、D3Q19、D3Q27分别是15、19、27个速度方向。D3Q15只有面心和角点缺了棱中点方向各向异性误差明显实际模拟中压力场容易变形现在用得很少。D3Q27把静止、面心、棱心、角点全带上精度上限最高但内存和计算量比D3Q19多出约四成而渗流模拟的工程精度通常用不到这个上限。D3Q19是折中——六个面心方向加十二个棱心方向正好覆盖三维空间的主轴和六个对角平面权重对称性好实现也不复杂。我做三维渗流基本默认D3Q19只有在孔隙结构极其复杂、压力梯度很高的场景才考虑D3Q27。三者的速度方向和权重有一个规律方向越短权重越大。D3Q19的权重是静止方向1/3、面心方向1/18、棱心方向1/36所有权重之和正好等于1。格子声速cs²在这个模型里恒等于1/3这是后面所有参数换算的基准。中国有句老话叫“差之毫厘谬以千里”权重表抄错一个数宏观方程就对不上了。我用的是下表这份速度模型方向构成稳定性内存开销适用场景D3Q15静止1 面心6 角点8最弱压力场易畸变最低教学演示不推荐工程D3Q19静止1 面心6 棱心12良好适中渗流、多孔介质、微流动D3Q27静止1 面心6 棱心12 角点8最优约为D3Q19的1.4倍高精度、强非平衡流动渗流模拟的雷诺数通常远小于1速度场变化平缓D3Q19的精度已经完全够用。真正影响结果的是后面的孔隙建模和边界处理而不是多两个速度方向。2.2 碰撞算子BGK够用什么时候换MRT三维LBM的碰撞算子最常用的是BGK也叫单松弛模型。它的核心思想很简单碰撞后的分布函数向平衡态松弛松弛速度由无量纲松弛时间tau决定。BGK的碰撞公式写成数组操作就一行f f - omega * (f - feq)其中omega 1 / taufeq是平衡态分布函数。这个式子对每个格子、每个方向独立计算天然适合向量化。BGK的问题是所有物理模态用同一个松弛速度在高雷诺数或高马赫数下容易激发非物理振荡。渗流场景下流速极低马赫数通常小于0.05BGK不会遇到这个问题。我一般只在tau被迫小于0.55、或者网格分辨率加密后仍出现振荡时才换MRT多松弛模型。MRT把分布函数投影到矩空间对不同矩用不同松弛系数数值稳定性明显更好但代价是代码量翻倍——每个格点要做两次19×19的矩阵变换三维数组的内存访问模式也变得不友好。对绝大多数孔隙渗流模拟BGK加体积力驱动已经足够稳定先把BGK跑通再根据实际表现决定要不要上MRT这是最务实的路径。2.3 松弛时间tau和网格分辨率决定收敛速度的两个旋钮tau是三维LBM里最敏感的参数。运动粘度nu与tau的关系是nu cs² * (tau - 0.5)把cs² 1/3代入得到nu (tau - 0.5) / 3。tau越接近0.5数值粘度越小物理上越接近无粘流动但计算稳定性越差tau太大数值耗散太大低速渗流的驱动力会被吃掉。我的经验区间是tau取0.6到0.8对应nu在0.033到0.1之间。渗流模拟中tau每增加0.1达到稳态需要的迭代步数可能多出两成。网格分辨率直接影响格子单位与物理单位的换算。三维渗流一般要求最小孔隙通道里至少保留2到3个格子来分辨近壁面边界层整体孔径最好能覆盖10个格子以上。格子分辨率不够时渗透率计算结果会明显偏离真实值而且这个误差是系统性的不是随机噪声。我通常的做法是先用32³或64³跑通流程再在最小孔径方向加密一倍做网格无关性验证两次结果差异小于2%才认定分辨率达标。3. 把渗流问题翻译成三维LBM孔隙建模、边界条件与渗透率标定3.1 多孔介质骨架生成随机概率法与四参数生长法三维渗流模拟的第一步是生成多孔介质骨架——用一个三维布尔数组标记每个格子是流体还是固体。最简单的生成方式是随机概率法每个格子以概率p被标记为固体1 - p就是目标孔隙度。这个方法的优点是代码只有几行缺点是孔隙连通性完全不受控孔隙度稍高一点就会出现大量孤立的固体团块或死端孔隙渗透率会严重低估真实多孔介质。我一般只拿它做算法调试不用于正式研究。import numpy as np from scipy import ndimage def generate_random_porous(nx, ny, nz, phi, seed42): 随机概率法生成多孔介质。 phi为目标孔隙度返回的solid数组中True代表固体格点。 rng np.random.default_rng(seed) solid rng.random((nx, ny, nz)) (1.0 - phi) # 保证出入口三层面全部为流体避免边界被堵 solid[0, :, :] False solid[-1, :, :] False return solid随机概率法的核心参数是孔隙度phi它直接决定固体格点的比例。代码最后两行把x方向的首尾两层强制设为流体是为了配合后面的体积力驱动边界否则入口被固体堵住渗流根本建立不起来。需要注意的是随机生成后要做连通性检查用ndimage.label统计流体区域的连通分量主渗流通道必须是连通的否则计算出的渗透率会是零或噪声。更接近真实岩心的做法是四参数生长法QSGS。它先以很小的核心概率p_c随机撒种子格点作为固体生长核心然后每个核心向6个相邻方向以生长概率p_i扩展重复多轮。这样生成的孔隙比随机法更连通渗透率随孔隙度的变化也更平滑。核心概率p_c通常取0.01到0.05p_i取0.05到0.3生长轮数决定孔隙的迂曲度。QSGS的参数不像随机法只有一个phi它需要根据目标孔隙度反推p_c和p_i组合第一次试跑时建议固定p_c扫描p_i并统计生成的孔隙度。3.2 边界条件取舍体积力驱动还是压力边界多孔介质渗流的边界条件有两条路压力边界驱动或体积力驱动。压力边界需要在入口和出口各施加恒定的密度差用Zou-He压力边界或非平衡外推实现代码复杂且在多孔介质这种非平坦边界上容易引入数值伪影。体积力驱动的做法是在每个流体格点上施加一个恒定的外力配合周期性边界自动形成压力梯度实现简单稳定性好被绝大多数孔隙尺度LBM研究采用。我做三维渗流一致推荐体积力驱动。周期性边界在三维LBM里用np.roll就能实现。np.roll会把数组沿指定方向循环平移从一边移出去的数据自动从另一边补回来天然满足周期条件。核心意思是每个方向上的分布函数迁移就是整层数组平移一格。这个操作效率很高但也埋了一个坑——np.roll会同时移动流体格点和固体格点的数据所以每步迁移后必须做反弹边界处理把进入固体的分布函数送回流体区域。固体边界的处理用半程反弹法。它的物理含义是分布函数碰到固体壁面后原路返回等效于壁面位于流体格点和固体格点正中间。实现时维护一个opp数组记录每个方向的反方向编号方向0对应自身方向1对应2方向3对应4以此类推。反弹操作在三维LBM里是唯一的难点也是翻车的高发区在第4章的代码里我会给出可用的实现。3.3 从速度场到渗透率稳态判据与达西定律渗流模拟跑到稳态后需要从速度场里提取渗透率。稳态判据最常用的是逐迭代步的全局速度变化率连续若干步内所有格点速度的最大变化小于一个阈值比如1e-9就认为到达稳态。更工程化的做法是同时监测出口截面的平均流速它收敛到常数且进出口流量差不大于1%时进入数据采样阶段。渗透率由达西定律反推。达西定律在体积力驱动的周期边界下写为k (Q / A) * rho * nu / F其中Q是过某截面的体积流量A是截面面积rho是流体密度nu是运动粘度F是施加的体积力强度。在LBM格子单位里rho约等于1nu由tau换算得到Q/A就是截面平均速度U_avg最后得到k_lu U_avg * nu / F这里算出的渗透率是格子单位换成物理渗透率要乘上格子尺寸dx的平方。1达西约等于0.987e-12平方米很多初学者算出来的渗透率数值比岩心实测值大好几个数量级原因就是忘了乘dx²。我做这类计算时习惯在代码里保留dx为显式参数从格子单位到物理单位的换算放最后统一处理。4. 用Python跑通三维LBM渗流可复现的D3Q19最小实现4.1 数据结构与平衡态分布先搭骨架三维LBM的分布函数数组我习惯做成f[..., n]的形状n是最后一个维度大小为19。这种内存布局的好处是np.roll沿空间轴平移时分布函数的各个分量连续排列cache命中率比方向在前的布局高不少对大网格的性能差别非常明显。固体标记数组solid做成布尔数组True代表固体计算时用float类型的mask做掩码。import numpy as np # D3Q19 速度方向0静止1-6面心7-18棱心 c np.array([ [ 0, 0, 0], # 0 [ 1, 0, 0], [-1, 0, 0], [ 0, 1, 0], [ 0,-1, 0], [ 0, 0, 1], [ 0, 0,-1], # 1-6 [ 1, 1, 0], [-1,-1, 0], [ 1,-1, 0], [-1, 1, 0], # 7-10 [ 1, 0, 1], [-1, 0,-1], [ 1, 0,-1], [-1, 0, 1], # 11-14 [ 0, 1, 1], [ 0,-1,-1], [ 0, 1,-1], [ 0,-1, 1], # 15-18 ], dtypenp.int64) w np.array([1/3] [1/18]*6 [1/36]*12, dtypenp.float64) opp np.array([0, 2, 1, 4, 3, 6, 5, 9, 8, 7, 10, 13, 12, 11, 14, 17, 16, 15, 18]) cs2 1.0 / 3.0这段代码定义了D3Q19的全部19个速度方向、对应的平衡态权重以及反弹边界需要的反方向索引。opp数组的构造逻辑是每个方向的反方向方向1向右反方向2向左方向7是(1,1,0)反方向8是(-1,-1,0)一一对应。这个数组一旦写错反弹边界会把分布函数送回错误的格点模拟结果立刻变成垃圾。我每次新建项目都会用np.allclose核对opp[c]恰好是-c这比肉眼检查可靠得多。平衡态分布函数是BGK碰撞的基础在低马赫数下写成def equilibrium(rho, ux, uy, uz): D3Q19的平衡态分布函数。 cs2 1/3公式为 f_eq w * rho * (1 3*(c·u) 4.5*(c·u)^2 - 1.5*u^2) usqr (ux**2 uy**2 uz**2).astype(np.float64) cu (ux[..., None] * c[:, 0] uy[..., None] * c[:, 1] uz[..., None] * c[:, 2]) feq w * rho[..., None] * (1.0 3.0 * cu 4.5 * cu**2 - 1.5 * usqr[..., None]) return feq平衡态公式里cu是每个格点上宏观速度在各自速度方向上的投影rho[..., None]把密度扩展到19个方向与权重w逐方向相乘。公式中的第二项对应压力项第三项和第四项对应动量通量缺了任何一项都会导致宏观方程不完整。常见错误是在cu**2那项把系数4.5写成1.5结果就是压力场怎么迭代都对不上解析解。4.2 碰撞、迁移与反弹边界主循环怎么组装主循环每步分为三段碰撞、迁移、反弹。碰撞阶段用BGK更新全部格点的分布函数迁移阶段用np.roll把分布平移到相邻格点反弹阶段把迁移进固体的分布送回流体。顺序不能颠倒这是半程反弹的标准流程。def collision(f, rho, ux, uy, uz, omega): BGK碰撞返回碰撞后的分布函数。 feq equilibrium(rho, ux, uy, uz) return f - omega * (f - feq) def streaming(f): 周期性迁移按速度方向平移分布函数。 f_new np.empty_like(f) for i in range(19): ci c[i] f_new[..., i] np.roll(f[..., i], shiftci, axis(0, 1, 2)) return f_new def bounce_back(f, solid): 半程反弹固体格点内迁移进来的分布函数按反方向送回流体。 solid_f solid.astype(np.float64) f[..., :] * (1.0 - solid_f)[..., None] # 先清空固体内的分布 for i in range(19): j opp[i] # 找出哪些格点的i方向邻居是固体把solid按反方向平移一格 src_mask np.roll(solid_f, shift-c[i], axis(0, 1, 2)) # 固体格点里的i方向分布弹回源格点的j方向 incoming f[..., i] * solid_f f[..., j] np.roll(incoming, shift-c[i], axis(0, 1, 2)) f[..., i] - incoming return f迁移函数里np.roll的shift传入三元组分别对应x、y、z三个轴的平移量一次调用完成三维平移。反弹函数的逻辑分三步先把固体格点里的所有分布函数清零再对每个方向i找出“目标邻居为固体”的流体源格点把已经迁移到固体里的数值取出来原路平移回源格点并加到反方向分量上最后从源方向里减掉。这样保证每步结束后固体格点里没有残留分布流量守恒不会被破坏。这个反弹实现有一个近似同一格点在同一时间步收到来自多个方向的反弹时会直接累加严格意义上应该在累加前固定反弹源值。但在渗流这种低马赫数场景误差远小于0.1%工程上可以接受。如果做高精度定量研究可以把反弹分拆到迁移阶段分别处理代价是代码可读性下降。4.3 渗流主程序从多孔介质到渗透率输出主程序把上面的函数组装起来加上体积力驱动项和宏观量提取。体积力用Guo力项它在低马赫数下比简单的“速度偏移法”更稳定能有效抑制棋盘振荡。主程序跑起来之后每100步记录一次截面平均速度速度变化小于设定阈值就停机。def guo_force(rho, ux, uy, uz, fx, omega): Guo体积力项适用于D3Q19。fx为x方向体积力。 cu (ux[..., None] * c[:, 0] uy[..., None] * c[:, 1] uz[..., None] * c[:, 2]) force np.zeros_like(fx, dtypenp.float64) force[..., None] fx[..., None] * c[:, 0] # 展成19个方向 s (1.0 - 0.5 * omega) * w * ( 3.0 * force 9.0 * cu * force - 3.0 * ux[..., None] * fx[..., None] ) return s def run_porous_lbm(solid, phi, tau, fx, steps20000, tol1e-8): nx, ny, nz solid.shape omega 1.0 / tau fluid (1.0 - solid.astype(np.float64))[..., None] f np.broadcast_to(w, (nx, ny, nz, 19)).copy() f * fluid # 固体格点初始为0 history [] for step in range(steps): rho f.sum(axis-1) 1e-12 ux (f[..., 1] - f[..., 2] f[..., 7] - f[..., 8] f[..., 11] - f[..., 12] f[..., 13] - f[..., 14]) / rho uy (f[..., 3] - f[..., 4] f[..., 7] - f[..., 8] f[..., 9] - f[..., 10] f[..., 15] - f[..., 16]) / rho uz (f[..., 5] - f[..., 6] f[..., 11] - f[..., 12] f[..., 13] - f[..., 14] f[..., 17] - f[..., 18]) / rho f collision(f, rho, ux, uy, uz, omega) f guo_force(rho, ux, uy, uz, fx, omega) f streaming(f) f bounce_back(f, solid.astype(np.float64)) if step % 100 0: ux_avg ux[(1 - solid).astype(bool)].mean() history.append(ux_avg) if len(history) 20 and abs(history[-1] - history[-20]) tol: break u_avg np.mean(ux[(1 - solid).astype(bool)]) nu cs2 * (tau - 0.5) k_lu u_avg * nu / fx return k_lu, ux, history这段代码把碰撞、外力、迁移、反弹整合进单个循环。宏观速度的计算只用到19个方向中带x、y、z分量的那些等价于用方向0以外的分布求和再除以密度。固体格点的rho被加了一个1e-12的保护项避免全零分布时出现除零。Guo力项的系数(1 - 0.5*omega)是体积力在BGK模型里的标准修正直接加力不加这个系数会导致渗透率偏低。稳态判断用的是每100步记录一次截面平均速度对比当前值与20个记录前的差值。溢出判断是第一个该加上的防线——任何一步出现np.isnan(ux)就要立即停止并打印当前step和tau而不是等它跑满20000步。渗透率的计算放最后k_lu乘上dx²才是物理渗透率主程序里把dx作为参数传入我一般写成k_phys k_lu * dx * dx。5. 三维LBM渗流排查四个高频翻车点与对应解法5.1 跑着跑着变成NaN松弛时间和力项首当其冲现象循环跑几千步后速度场和密度场突然出现NaN整个数组在几十步内全部污染。原因最常见的是tau取值过小tau接近0.5时数值粘度接近零任何微小扰动都会被放大。另一个高频原因是体积力fx设置过大导致局部格点的马赫数超过0.1非线性项失去数值约束。还有可能是Guo力项里的宏观速度u用了碰撞前的值力项与碰撞不同步形成正反馈振荡。解决先把tau升到0.7重新跑如果稳定再把tau逐步降到0.6。fx的初始值我用0.0001量级观察截面平均速度接近0.01量级后再按比例增大。力项差值的检查方法是在循环里单独输出固体边界附近的密度最大值如果出现超过1.1的格点说明反弹边界或力项附近速度过冲需要把fx减半再验证。5.2 速度场出棋盘振荡边界mask处理不干净现象三维速度切面图上出现红绿相间的棋盘状斑块尤其是在固体与流体交界附近速度方向交替反转。原因棋盘振荡是LBM的经典数值病态通常来自力项处理与碰撞的不一致。如果用了速度偏移法施加外力而平衡态又用了碰撞前的宏观速度压力信息会以锯齿方式传播。另一个隐蔽原因是反弹边界里固体mask没有清零干净固体格点残留分布被周期迁移带到流体区域形成持续的高频扰动。解决统一改用Guo力项并确认系数(1 - 0.5*omega)没有写错。检查bounce_back函数里9到12行确保进入固体的每个方向都被减去而不是只加了反方向。验证方法是在纯流体空场无固体下跑同体积力理论上速度场必须严格均匀且无振荡如果空场都振荡问题出在力项而不是边界。5.3 渗透率不随网格收敛边界层分辨率不够现象把网格从32³加密到64³再加密到96³渗透率一直在变没有收敛趋势。通常是小网格算出的渗透率偏大加密后明显下降。原因这是典型的近壁面边界层未解析。半程反弹边界等效于把壁面放在流体格子和固体格子中点当孔隙通道只有三四个格子宽时反弹位置与真实壁面的偏差占孔隙宽度的比例过大相当于把孔隙人为加宽了。孔径越小这个虚增的通道越宽渗透率就越高。工程上把这个现象归因于格子Knudsen数效应——格子分辨率不够时边界层内的速度梯度根本铺不开。解决做法是保证最小孔径方向至少有10个格子且边界层内至少有2到3个格子。我做网格无关性验证时每个方向独立加密一倍比较渗透率变化。如果加密后渗透率下降超过2%就继续加密已经到100³以上变化仍然明显说明多孔介质生成参数本身有问题——孔隙度太低导致连通通道太少建议重生成介质而不是盲目加密。另外tau取0.8时数值耗散会把边界层抹平必要时把tau降到0.6再看收敛性。5.4 算得实在太慢内存布局和迁移写法的锅现象64³网格跑2万步要几个小时比同样规模的二维模拟慢了几十倍完全没法做参数扫描。原因纯Python循环的逐方向迁移是性能瓶颈。如果分布数组用f[19, nx, ny, nz]布局np.roll时要跨大跨度内存访问cache命中率极低。另一个常见坑是在碰撞函数里用Python循环遍历19个方向逐方向计算feq这种写法在NumPy里虽然能跑但会反复创建中间数组内存带宽被白白浪费。解决分布数组改用f[nx, ny, nz, 19]布局碰撞和迁移都用向量化表达式一次处理全数组。更进一步的方案是检查np.roll的调用次数——每步19次np.roll是必须的但np.roll返回新数组会复制数据可以改成预分配f_new逐方向写入。如果64³以上的网格还要做批量参数扫描建议直接切换到Taichi或CUDA实现同样的算法在GPU上通常能带来30到80倍加速纯Python优化的天花板有限。6. 进阶验证与提速可视化、二三维联动和GPU改造6.1 三维速度场可视化确认渗流通道是否连通三维LBM跑的到底对不对光看渗透率数值不够还要看速度场形态。我一般用matplotlib画出三维速度等值面或者在主渗流方向截几个切片把速度大小叠加到孔隙结构上。判断标准很简单主通路上的速度必须连贯固体周边速度趋近于零死端孔隙内速度显著低于主通道。如果等值面被割裂成互不相连的碎块多半是孔隙连通性检查没做生成的骨架本身就是断的。可视化这一步不花多少时间但能拦住后面一大堆无效计算。6.2 二三维联动标定用2D小模型定参数用3D出结果三维LBM的调参成本远高于二维一个64³算例就要几分钟。我的习惯是做二三维联动标定先用二维D2Q9模型在同一套多孔介质生成参数下快速扫描tau和fx找到稳定区间后再把参数搬到三维跑正式算例。这样做的前提是二维模型只用于定参数趋势不直接作为渗透率结果。二维与三维的渗透率天然有差异——二维把孔隙近似成无限长柱体迂曲度分布跟三维不一样所以最终渗透率必须以三维为准。二维的作用是省时间、快速排除明显不稳定的参数组合全钒液流电池多孔电极这类工程对象的渗流模拟我也沿用这套标定流程。6.3 GPU提速的切入点先动迁移再动碰撞当三维网格超过128³时纯NumPy已经撑不住。GPU改造不必一步到位我建议先改迁移再改碰撞。迁移是访问模式最规整的阶段——每个格点向19个方向搬运数据天然并行碰撞也同理每个格点独立计算。用CuPy时np.roll可以几乎不改动地换成cp.roll分布函数数组换成Cupy数组等待编译的时间比代码改写的时间还长。Taichi的表现更稳速度方向和权重用ti.static声明迁移循环可以写成逐格点并行不受NumPy数组shape限制。实际经验是先确认反弹边界的mask逻辑在GPU上完全一致后再考虑优化力项和平衡态函数的寄存器占用。每次缩小半分冤枉路整套流程跑通后三维渗流模拟才能真正摆脱“跑几步就等半天”的尴尬。如果你手头的需求是页岩渗流、土壤污染物迁移或多孔电极内电解液流动三维LBM这套路径都适用。我自己的教训是多孔介质生成阶段多花一小时做连通性检查比在LBM主循环里排查一整天要划算得多。参数调不明白时先回到最小网格纯流体算例把碰撞、迁移、反弹逐段验证再叠加多孔介质。希望这篇笔记能让你少走几段弯路把三维LBM渗流的第一个算例顺利跑出来。本文还有配套的精品资源点击获取
简介这是一份基于格子玻尔兹曼方法LBM的三维渗流模拟MATLAB源码主要面向流体力学、多孔介质渗流方向的工程师与科研人员适用于地下水流动、石油开采、土壤污染扩散等典型渗流场景。该方法通过离散玻尔兹曼方程模拟微观粒子运动进而恢复宏观流动行为对复杂边界和非均匀介质有天然优势三维算法相比二维模拟能呈现更全面的流场信息更贴近工程实际。压缩包为zip格式共1个文件为m类型脚本整体仅3KB属于轻量级核心代码方便直接阅读算法流程与参数配置。目前已有601人学习下载。借助该源码使用者可自定义流体参数、网格尺寸与边界条件执行初始化分布函数、碰撞、迁移、边界处理等标准步骤最终导出流速、压力分布的可视化流场图并计算无量纲渗流系数为渗流特性分析提供重要数据支撑。代码结构简洁既适合作为三维LBM渗流教学演示也可在此基础上做科研二次开发。1. 三维LBM渗流孔隙尺度的流动模拟入门到底难在哪“LBM三维”和“LBM渗流”连在一起指向的是一类很具体的需求用格子玻尔兹曼方法Lattice Boltzmann MethodLBM在孔隙尺度上把多孔介质中的渗流过程算出来。跟达西尺度的黑匣子模型不同三维LBM能还原每个孔隙里的速度场和压力场天然适配岩芯、页岩、土壤、多孔电极这类非均质结构。对刚入门的人它是一套可以完全用NumPy搭起来的数值框架对从二维转向三维的工程师最大的成本不在算法而在模型选型、边界处理和稳定性调参这三件事。这篇文章按我实际跑通的经验把D3Q19模型、多孔介质生成、反弹边界到渗透率输出这条链路拆开讲代码可以直接抄。2. 三维LBM选型D3Q19、松弛时间与稳定性边界2.1 三维速度模型怎么选D3Q15、D3Q19、D3Q27的取舍三维LBM比二维多出来的第一个选择就是离散速度模型。D3Q15、D3Q19、D3Q27分别是15、19、27个速度方向。D3Q15只有面心和角点缺了棱中点方向各向异性误差明显实际模拟中压力场容易变形现在用得很少。D3Q27把静止、面心、棱心、角点全带上精度上限最高但内存和计算量比D3Q19多出约四成而渗流模拟的工程精度通常用不到这个上限。D3Q19是折中——六个面心方向加十二个棱心方向正好覆盖三维空间的主轴和六个对角平面权重对称性好实现也不复杂。我做三维渗流基本默认D3Q19只有在孔隙结构极其复杂、压力梯度很高的场景才考虑D3Q27。三者的速度方向和权重有一个规律方向越短权重越大。D3Q19的权重是静止方向1/3、面心方向1/18、棱心方向1/36所有权重之和正好等于1。格子声速cs²在这个模型里恒等于1/3这是后面所有参数换算的基准。中国有句老话叫“差之毫厘谬以千里”权重表抄错一个数宏观方程就对不上了。我用的是下表这份速度模型方向构成稳定性内存开销适用场景D3Q15静止1 面心6 角点8最弱压力场易畸变最低教学演示不推荐工程D3Q19静止1 面心6 棱心12良好适中渗流、多孔介质、微流动D3Q27静止1 面心6 棱心12 角点8最优约为D3Q19的1.4倍高精度、强非平衡流动渗流模拟的雷诺数通常远小于1速度场变化平缓D3Q19的精度已经完全够用。真正影响结果的是后面的孔隙建模和边界处理而不是多两个速度方向。2.2 碰撞算子BGK够用什么时候换MRT三维LBM的碰撞算子最常用的是BGK也叫单松弛模型。它的核心思想很简单碰撞后的分布函数向平衡态松弛松弛速度由无量纲松弛时间tau决定。BGK的碰撞公式写成数组操作就一行f f - omega * (f - feq)其中omega 1 / taufeq是平衡态分布函数。这个式子对每个格子、每个方向独立计算天然适合向量化。BGK的问题是所有物理模态用同一个松弛速度在高雷诺数或高马赫数下容易激发非物理振荡。渗流场景下流速极低马赫数通常小于0.05BGK不会遇到这个问题。我一般只在tau被迫小于0.55、或者网格分辨率加密后仍出现振荡时才换MRT多松弛模型。MRT把分布函数投影到矩空间对不同矩用不同松弛系数数值稳定性明显更好但代价是代码量翻倍——每个格点要做两次19×19的矩阵变换三维数组的内存访问模式也变得不友好。对绝大多数孔隙渗流模拟BGK加体积力驱动已经足够稳定先把BGK跑通再根据实际表现决定要不要上MRT这是最务实的路径。2.3 松弛时间tau和网格分辨率决定收敛速度的两个旋钮tau是三维LBM里最敏感的参数。运动粘度nu与tau的关系是nu cs² * (tau - 0.5)把cs² 1/3代入得到nu (tau - 0.5) / 3。tau越接近0.5数值粘度越小物理上越接近无粘流动但计算稳定性越差tau太大数值耗散太大低速渗流的驱动力会被吃掉。我的经验区间是tau取0.6到0.8对应nu在0.033到0.1之间。渗流模拟中tau每增加0.1达到稳态需要的迭代步数可能多出两成。网格分辨率直接影响格子单位与物理单位的换算。三维渗流一般要求最小孔隙通道里至少保留2到3个格子来分辨近壁面边界层整体孔径最好能覆盖10个格子以上。格子分辨率不够时渗透率计算结果会明显偏离真实值而且这个误差是系统性的不是随机噪声。我通常的做法是先用32³或64³跑通流程再在最小孔径方向加密一倍做网格无关性验证两次结果差异小于2%才认定分辨率达标。3. 把渗流问题翻译成三维LBM孔隙建模、边界条件与渗透率标定3.1 多孔介质骨架生成随机概率法与四参数生长法三维渗流模拟的第一步是生成多孔介质骨架——用一个三维布尔数组标记每个格子是流体还是固体。最简单的生成方式是随机概率法每个格子以概率p被标记为固体1 - p就是目标孔隙度。这个方法的优点是代码只有几行缺点是孔隙连通性完全不受控孔隙度稍高一点就会出现大量孤立的固体团块或死端孔隙渗透率会严重低估真实多孔介质。我一般只拿它做算法调试不用于正式研究。import numpy as np from scipy import ndimage def generate_random_porous(nx, ny, nz, phi, seed42): 随机概率法生成多孔介质。 phi为目标孔隙度返回的solid数组中True代表固体格点。 rng np.random.default_rng(seed) solid rng.random((nx, ny, nz)) (1.0 - phi) # 保证出入口三层面全部为流体避免边界被堵 solid[0, :, :] False solid[-1, :, :] False return solid随机概率法的核心参数是孔隙度phi它直接决定固体格点的比例。代码最后两行把x方向的首尾两层强制设为流体是为了配合后面的体积力驱动边界否则入口被固体堵住渗流根本建立不起来。需要注意的是随机生成后要做连通性检查用ndimage.label统计流体区域的连通分量主渗流通道必须是连通的否则计算出的渗透率会是零或噪声。更接近真实岩心的做法是四参数生长法QSGS。它先以很小的核心概率p_c随机撒种子格点作为固体生长核心然后每个核心向6个相邻方向以生长概率p_i扩展重复多轮。这样生成的孔隙比随机法更连通渗透率随孔隙度的变化也更平滑。核心概率p_c通常取0.01到0.05p_i取0.05到0.3生长轮数决定孔隙的迂曲度。QSGS的参数不像随机法只有一个phi它需要根据目标孔隙度反推p_c和p_i组合第一次试跑时建议固定p_c扫描p_i并统计生成的孔隙度。3.2 边界条件取舍体积力驱动还是压力边界多孔介质渗流的边界条件有两条路压力边界驱动或体积力驱动。压力边界需要在入口和出口各施加恒定的密度差用Zou-He压力边界或非平衡外推实现代码复杂且在多孔介质这种非平坦边界上容易引入数值伪影。体积力驱动的做法是在每个流体格点上施加一个恒定的外力配合周期性边界自动形成压力梯度实现简单稳定性好被绝大多数孔隙尺度LBM研究采用。我做三维渗流一致推荐体积力驱动。周期性边界在三维LBM里用np.roll就能实现。np.roll会把数组沿指定方向循环平移从一边移出去的数据自动从另一边补回来天然满足周期条件。核心意思是每个方向上的分布函数迁移就是整层数组平移一格。这个操作效率很高但也埋了一个坑——np.roll会同时移动流体格点和固体格点的数据所以每步迁移后必须做反弹边界处理把进入固体的分布函数送回流体区域。固体边界的处理用半程反弹法。它的物理含义是分布函数碰到固体壁面后原路返回等效于壁面位于流体格点和固体格点正中间。实现时维护一个opp数组记录每个方向的反方向编号方向0对应自身方向1对应2方向3对应4以此类推。反弹操作在三维LBM里是唯一的难点也是翻车的高发区在第4章的代码里我会给出可用的实现。3.3 从速度场到渗透率稳态判据与达西定律渗流模拟跑到稳态后需要从速度场里提取渗透率。稳态判据最常用的是逐迭代步的全局速度变化率连续若干步内所有格点速度的最大变化小于一个阈值比如1e-9就认为到达稳态。更工程化的做法是同时监测出口截面的平均流速它收敛到常数且进出口流量差不大于1%时进入数据采样阶段。渗透率由达西定律反推。达西定律在体积力驱动的周期边界下写为k (Q / A) * rho * nu / F其中Q是过某截面的体积流量A是截面面积rho是流体密度nu是运动粘度F是施加的体积力强度。在LBM格子单位里rho约等于1nu由tau换算得到Q/A就是截面平均速度U_avg最后得到k_lu U_avg * nu / F这里算出的渗透率是格子单位换成物理渗透率要乘上格子尺寸dx的平方。1达西约等于0.987e-12平方米很多初学者算出来的渗透率数值比岩心实测值大好几个数量级原因就是忘了乘dx²。我做这类计算时习惯在代码里保留dx为显式参数从格子单位到物理单位的换算放最后统一处理。4. 用Python跑通三维LBM渗流可复现的D3Q19最小实现4.1 数据结构与平衡态分布先搭骨架三维LBM的分布函数数组我习惯做成f[..., n]的形状n是最后一个维度大小为19。这种内存布局的好处是np.roll沿空间轴平移时分布函数的各个分量连续排列cache命中率比方向在前的布局高不少对大网格的性能差别非常明显。固体标记数组solid做成布尔数组True代表固体计算时用float类型的mask做掩码。import numpy as np # D3Q19 速度方向0静止1-6面心7-18棱心 c np.array([ [ 0, 0, 0], # 0 [ 1, 0, 0], [-1, 0, 0], [ 0, 1, 0], [ 0,-1, 0], [ 0, 0, 1], [ 0, 0,-1], # 1-6 [ 1, 1, 0], [-1,-1, 0], [ 1,-1, 0], [-1, 1, 0], # 7-10 [ 1, 0, 1], [-1, 0,-1], [ 1, 0,-1], [-1, 0, 1], # 11-14 [ 0, 1, 1], [ 0,-1,-1], [ 0, 1,-1], [ 0,-1, 1], # 15-18 ], dtypenp.int64) w np.array([1/3] [1/18]*6 [1/36]*12, dtypenp.float64) opp np.array([0, 2, 1, 4, 3, 6, 5, 9, 8, 7, 10, 13, 12, 11, 14, 17, 16, 15, 18]) cs2 1.0 / 3.0这段代码定义了D3Q19的全部19个速度方向、对应的平衡态权重以及反弹边界需要的反方向索引。opp数组的构造逻辑是每个方向的反方向方向1向右反方向2向左方向7是(1,1,0)反方向8是(-1,-1,0)一一对应。这个数组一旦写错反弹边界会把分布函数送回错误的格点模拟结果立刻变成垃圾。我每次新建项目都会用np.allclose核对opp[c]恰好是-c这比肉眼检查可靠得多。平衡态分布函数是BGK碰撞的基础在低马赫数下写成def equilibrium(rho, ux, uy, uz): D3Q19的平衡态分布函数。 cs2 1/3公式为 f_eq w * rho * (1 3*(c·u) 4.5*(c·u)^2 - 1.5*u^2) usqr (ux**2 uy**2 uz**2).astype(np.float64) cu (ux[..., None] * c[:, 0] uy[..., None] * c[:, 1] uz[..., None] * c[:, 2]) feq w * rho[..., None] * (1.0 3.0 * cu 4.5 * cu**2 - 1.5 * usqr[..., None]) return feq平衡态公式里cu是每个格点上宏观速度在各自速度方向上的投影rho[..., None]把密度扩展到19个方向与权重w逐方向相乘。公式中的第二项对应压力项第三项和第四项对应动量通量缺了任何一项都会导致宏观方程不完整。常见错误是在cu**2那项把系数4.5写成1.5结果就是压力场怎么迭代都对不上解析解。4.2 碰撞、迁移与反弹边界主循环怎么组装主循环每步分为三段碰撞、迁移、反弹。碰撞阶段用BGK更新全部格点的分布函数迁移阶段用np.roll把分布平移到相邻格点反弹阶段把迁移进固体的分布送回流体。顺序不能颠倒这是半程反弹的标准流程。def collision(f, rho, ux, uy, uz, omega): BGK碰撞返回碰撞后的分布函数。 feq equilibrium(rho, ux, uy, uz) return f - omega * (f - feq) def streaming(f): 周期性迁移按速度方向平移分布函数。 f_new np.empty_like(f) for i in range(19): ci c[i] f_new[..., i] np.roll(f[..., i], shiftci, axis(0, 1, 2)) return f_new def bounce_back(f, solid): 半程反弹固体格点内迁移进来的分布函数按反方向送回流体。 solid_f solid.astype(np.float64) f[..., :] * (1.0 - solid_f)[..., None] # 先清空固体内的分布 for i in range(19): j opp[i] # 找出哪些格点的i方向邻居是固体把solid按反方向平移一格 src_mask np.roll(solid_f, shift-c[i], axis(0, 1, 2)) # 固体格点里的i方向分布弹回源格点的j方向 incoming f[..., i] * solid_f f[..., j] np.roll(incoming, shift-c[i], axis(0, 1, 2)) f[..., i] - incoming return f迁移函数里np.roll的shift传入三元组分别对应x、y、z三个轴的平移量一次调用完成三维平移。反弹函数的逻辑分三步先把固体格点里的所有分布函数清零再对每个方向i找出“目标邻居为固体”的流体源格点把已经迁移到固体里的数值取出来原路平移回源格点并加到反方向分量上最后从源方向里减掉。这样保证每步结束后固体格点里没有残留分布流量守恒不会被破坏。这个反弹实现有一个近似同一格点在同一时间步收到来自多个方向的反弹时会直接累加严格意义上应该在累加前固定反弹源值。但在渗流这种低马赫数场景误差远小于0.1%工程上可以接受。如果做高精度定量研究可以把反弹分拆到迁移阶段分别处理代价是代码可读性下降。4.3 渗流主程序从多孔介质到渗透率输出主程序把上面的函数组装起来加上体积力驱动项和宏观量提取。体积力用Guo力项它在低马赫数下比简单的“速度偏移法”更稳定能有效抑制棋盘振荡。主程序跑起来之后每100步记录一次截面平均速度速度变化小于设定阈值就停机。def guo_force(rho, ux, uy, uz, fx, omega): Guo体积力项适用于D3Q19。fx为x方向体积力。 cu (ux[..., None] * c[:, 0] uy[..., None] * c[:, 1] uz[..., None] * c[:, 2]) force np.zeros_like(fx, dtypenp.float64) force[..., None] fx[..., None] * c[:, 0] # 展成19个方向 s (1.0 - 0.5 * omega) * w * ( 3.0 * force 9.0 * cu * force - 3.0 * ux[..., None] * fx[..., None] ) return s def run_porous_lbm(solid, phi, tau, fx, steps20000, tol1e-8): nx, ny, nz solid.shape omega 1.0 / tau fluid (1.0 - solid.astype(np.float64))[..., None] f np.broadcast_to(w, (nx, ny, nz, 19)).copy() f * fluid # 固体格点初始为0 history [] for step in range(steps): rho f.sum(axis-1) 1e-12 ux (f[..., 1] - f[..., 2] f[..., 7] - f[..., 8] f[..., 11] - f[..., 12] f[..., 13] - f[..., 14]) / rho uy (f[..., 3] - f[..., 4] f[..., 7] - f[..., 8] f[..., 9] - f[..., 10] f[..., 15] - f[..., 16]) / rho uz (f[..., 5] - f[..., 6] f[..., 11] - f[..., 12] f[..., 13] - f[..., 14] f[..., 17] - f[..., 18]) / rho f collision(f, rho, ux, uy, uz, omega) f guo_force(rho, ux, uy, uz, fx, omega) f streaming(f) f bounce_back(f, solid.astype(np.float64)) if step % 100 0: ux_avg ux[(1 - solid).astype(bool)].mean() history.append(ux_avg) if len(history) 20 and abs(history[-1] - history[-20]) tol: break u_avg np.mean(ux[(1 - solid).astype(bool)]) nu cs2 * (tau - 0.5) k_lu u_avg * nu / fx return k_lu, ux, history这段代码把碰撞、外力、迁移、反弹整合进单个循环。宏观速度的计算只用到19个方向中带x、y、z分量的那些等价于用方向0以外的分布求和再除以密度。固体格点的rho被加了一个1e-12的保护项避免全零分布时出现除零。Guo力项的系数(1 - 0.5*omega)是体积力在BGK模型里的标准修正直接加力不加这个系数会导致渗透率偏低。稳态判断用的是每100步记录一次截面平均速度对比当前值与20个记录前的差值。溢出判断是第一个该加上的防线——任何一步出现np.isnan(ux)就要立即停止并打印当前step和tau而不是等它跑满20000步。渗透率的计算放最后k_lu乘上dx²才是物理渗透率主程序里把dx作为参数传入我一般写成k_phys k_lu * dx * dx。5. 三维LBM渗流排查四个高频翻车点与对应解法5.1 跑着跑着变成NaN松弛时间和力项首当其冲现象循环跑几千步后速度场和密度场突然出现NaN整个数组在几十步内全部污染。原因最常见的是tau取值过小tau接近0.5时数值粘度接近零任何微小扰动都会被放大。另一个高频原因是体积力fx设置过大导致局部格点的马赫数超过0.1非线性项失去数值约束。还有可能是Guo力项里的宏观速度u用了碰撞前的值力项与碰撞不同步形成正反馈振荡。解决先把tau升到0.7重新跑如果稳定再把tau逐步降到0.6。fx的初始值我用0.0001量级观察截面平均速度接近0.01量级后再按比例增大。力项差值的检查方法是在循环里单独输出固体边界附近的密度最大值如果出现超过1.1的格点说明反弹边界或力项附近速度过冲需要把fx减半再验证。5.2 速度场出棋盘振荡边界mask处理不干净现象三维速度切面图上出现红绿相间的棋盘状斑块尤其是在固体与流体交界附近速度方向交替反转。原因棋盘振荡是LBM的经典数值病态通常来自力项处理与碰撞的不一致。如果用了速度偏移法施加外力而平衡态又用了碰撞前的宏观速度压力信息会以锯齿方式传播。另一个隐蔽原因是反弹边界里固体mask没有清零干净固体格点残留分布被周期迁移带到流体区域形成持续的高频扰动。解决统一改用Guo力项并确认系数(1 - 0.5*omega)没有写错。检查bounce_back函数里9到12行确保进入固体的每个方向都被减去而不是只加了反方向。验证方法是在纯流体空场无固体下跑同体积力理论上速度场必须严格均匀且无振荡如果空场都振荡问题出在力项而不是边界。5.3 渗透率不随网格收敛边界层分辨率不够现象把网格从32³加密到64³再加密到96³渗透率一直在变没有收敛趋势。通常是小网格算出的渗透率偏大加密后明显下降。原因这是典型的近壁面边界层未解析。半程反弹边界等效于把壁面放在流体格子和固体格子中点当孔隙通道只有三四个格子宽时反弹位置与真实壁面的偏差占孔隙宽度的比例过大相当于把孔隙人为加宽了。孔径越小这个虚增的通道越宽渗透率就越高。工程上把这个现象归因于格子Knudsen数效应——格子分辨率不够时边界层内的速度梯度根本铺不开。解决做法是保证最小孔径方向至少有10个格子且边界层内至少有2到3个格子。我做网格无关性验证时每个方向独立加密一倍比较渗透率变化。如果加密后渗透率下降超过2%就继续加密已经到100³以上变化仍然明显说明多孔介质生成参数本身有问题——孔隙度太低导致连通通道太少建议重生成介质而不是盲目加密。另外tau取0.8时数值耗散会把边界层抹平必要时把tau降到0.6再看收敛性。5.4 算得实在太慢内存布局和迁移写法的锅现象64³网格跑2万步要几个小时比同样规模的二维模拟慢了几十倍完全没法做参数扫描。原因纯Python循环的逐方向迁移是性能瓶颈。如果分布数组用f[19, nx, ny, nz]布局np.roll时要跨大跨度内存访问cache命中率极低。另一个常见坑是在碰撞函数里用Python循环遍历19个方向逐方向计算feq这种写法在NumPy里虽然能跑但会反复创建中间数组内存带宽被白白浪费。解决分布数组改用f[nx, ny, nz, 19]布局碰撞和迁移都用向量化表达式一次处理全数组。更进一步的方案是检查np.roll的调用次数——每步19次np.roll是必须的但np.roll返回新数组会复制数据可以改成预分配f_new逐方向写入。如果64³以上的网格还要做批量参数扫描建议直接切换到Taichi或CUDA实现同样的算法在GPU上通常能带来30到80倍加速纯Python优化的天花板有限。6. 进阶验证与提速可视化、二三维联动和GPU改造6.1 三维速度场可视化确认渗流通道是否连通三维LBM跑的到底对不对光看渗透率数值不够还要看速度场形态。我一般用matplotlib画出三维速度等值面或者在主渗流方向截几个切片把速度大小叠加到孔隙结构上。判断标准很简单主通路上的速度必须连贯固体周边速度趋近于零死端孔隙内速度显著低于主通道。如果等值面被割裂成互不相连的碎块多半是孔隙连通性检查没做生成的骨架本身就是断的。可视化这一步不花多少时间但能拦住后面一大堆无效计算。6.2 二三维联动标定用2D小模型定参数用3D出结果三维LBM的调参成本远高于二维一个64³算例就要几分钟。我的习惯是做二三维联动标定先用二维D2Q9模型在同一套多孔介质生成参数下快速扫描tau和fx找到稳定区间后再把参数搬到三维跑正式算例。这样做的前提是二维模型只用于定参数趋势不直接作为渗透率结果。二维与三维的渗透率天然有差异——二维把孔隙近似成无限长柱体迂曲度分布跟三维不一样所以最终渗透率必须以三维为准。二维的作用是省时间、快速排除明显不稳定的参数组合全钒液流电池多孔电极这类工程对象的渗流模拟我也沿用这套标定流程。6.3 GPU提速的切入点先动迁移再动碰撞当三维网格超过128³时纯NumPy已经撑不住。GPU改造不必一步到位我建议先改迁移再改碰撞。迁移是访问模式最规整的阶段——每个格点向19个方向搬运数据天然并行碰撞也同理每个格点独立计算。用CuPy时np.roll可以几乎不改动地换成cp.roll分布函数数组换成Cupy数组等待编译的时间比代码改写的时间还长。Taichi的表现更稳速度方向和权重用ti.static声明迁移循环可以写成逐格点并行不受NumPy数组shape限制。实际经验是先确认反弹边界的mask逻辑在GPU上完全一致后再考虑优化力项和平衡态函数的寄存器占用。每次缩小半分冤枉路整套流程跑通后三维渗流模拟才能真正摆脱“跑几步就等半天”的尴尬。如果你手头的需求是页岩渗流、土壤污染物迁移或多孔电极内电解液流动三维LBM这套路径都适用。我自己的教训是多孔介质生成阶段多花一小时做连通性检查比在LBM主循环里排查一整天要划算得多。参数调不明白时先回到最小网格纯流体算例把碰撞、迁移、反弹逐段验证再叠加多孔介质。希望这篇笔记能让你少走几段弯路把三维LBM渗流的第一个算例顺利跑出来。本文还有配套的精品资源点击获取
