简介一份基于Shan-Chen模型的3D多孔介质流动LBM仿真脚本面向流体力学、渗流力学及地质工程领域的研究人员与学生。Shan-Chen模型通过引入势能函数描述流体间相互作用能够处理多组分/多相流等复杂情形脚本以MATLAB实现可帮助读者直观获取多孔介质内速度场、渗透率与压力梯度分布适用于地下水流动、污染物扩散及储层渗流等场景的初步模拟与教学演示。包体共1个文件为.m脚本压缩包仅2KB代码精简便于阅读、修改与扩展。当前已有273人学习下载。借助该脚本和配套说明读者可深入理解Shan-Chen模型在三维网格上的离散与迭代实现掌握LBM多孔介质建模的基本流程并可作为OpenLB等其他LBM平台的对照参考。1. 为什么要把 Shan-Chen 模型搬进 3D 多孔介质从 LBM 到伪势力这个标题浓缩的是一个典型计算流体课题用格子 Boltzmann 方法LBM里的 Shan-Chen 伪势模型在 MATLAB 中实现对三维多孔介质内两相流动的模拟。Shan-Chen 模型不显式追踪界面而是用伪势在密度梯度上自动生成相分离和表面张力多孔介质中复杂的固体骨架又能用 LBM 的半程反弹格式相对自然地处理。三者合在一起能回答一个工程问题给定一块 3D 孔隙结构两相驱替的前缘怎么推进、剩余液滴卡在哪里、小孔道会不会被堵。适合油气渗流、燃料电池气体扩散层、岩心驱替实验的前期筛选。难点通常集中在三处3D 离散后的伪势力怎么接进演化方程、G 和 τ 这类参数怎么定、以及 MATLAB 版本怎样才能写得既正确又不至于慢到跑不动。2. D3Q19 离散、伪势力计算与 Shan-Chen 力模型的接入2.1 为什么多孔介质模拟默认选 D3Q19LBM 的三维速度集不止一种。D3Q19 表示三维空间里 19 个离散速度方向包含 1 个静止方向、6 个面心方向和 12 个棱心方向D3Q27 还会多出 8 个体心方向。教学代码和多数开源实现选 D3Q19主要原因是内存少约 30%碰撞和迁移的计算量明显更低而宏观方程不受影响。代价是对界面曲率的各向异性稍微敏感但 Shan-Chen 模型在多孔介质里的主要误差来源往往是伪速度spurious current不是速度集本身。所以先用 D3Q19 跑通流程确认两相分离和驱替曲线正常再考虑升级到 D3Q27。D3Q19 的常数必须和方向保持一一对应表 1 是方向分组和权重方向类别方向数速度方向示例权重 w_i静止1(0,0,0)1/3面心6(±1,0,0)、(0,±1,0)、(0,0,±1)1/18棱心12(±1,±1,0)、(±1,0,±1)、(0,±1,±1)1/36这类常数建议写成独立函数避免在演化代码里手滑打错。常见做法是function [c,w] d3q19() c [ 0 0 0; 1 0 0; -1 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 -1; 1 1 0; 1 -1 0; -1 1 0; -1 -1 0; 1 0 1; 1 0 -1; -1 0 1; -1 0 -1; 0 1 1; 0 1 -1; 0 -1 1; 0 -1 -1 ]; w [1/3; repmat(1/18,6,1); repmat(1/36,12,1)]; end这里把 c 写成 19×3 矩阵w 写成 19×1 列向量后续计算伪势力和平衡态分布时直接用c(i,:)和w(i)索引即可。注意第 1 行必须是静止速度否则权重对应关系会整体错位宏观量也会出错。2.2 伪势力流体-流体与流体-固体两个循环Shan-Chen 模型的核心是伪势力。流体-流体伪势力的标准写法是F_f(x) -G ψ(ρ(x)) ∑ w_i ψ(ρ(xe_i)) e_i其中 ψ 称为伪势函数G 控制两相之间的作用强度。G 为正时界面两侧密度互相排斥两相自动分离G 为负时表现为吸引。多孔介质两相驱替里几乎都用正 G。ψ(ρ) 可以取 ρ也可以取 ρ0(1-exp(-ρ/ρ0))。线性伪势实现简单但高密度区力增长太快密度稍高就容易数值震荡指数形式在大密度下趋向饱和稳定性明显更好。教学实现里指数形式更常见。固体骨架的润湿性通过流体-固体力接入常见形式是F_s(x) -G_s ψ(ρ(x)) ∑ w_i s(xe_i) e_i其中 s(x) 是固体指示函数固体节点取 1流体节点取 0。G_s 的正负和大小共同决定接触角取负时流体被固体吸引表现为亲液取正时表现为疏液。这个公式和流体-流体伪势力共用同一个循环只是在求和时额外乘一个邻居固体掩码。实际代码可以写成function F shanchen_force(rho, solid, G, Gs) % 计算 Shan-Chen 伪势力返回 Nx*Ny*Nz*3 的力场 % rho: 密度场solid: 逻辑数组固体为 true [c, w] d3q19(); psi 1 - exp(-rho); % 指数型伪势rho01 F zeros([size(rho), 3]); for i 2:19 rhoN circshift(psi, -c(i,:)); % 邻居密度 sN circshift(solid, -c(i,:)); % 邻居固体标记 fcom G * psi .* rhoN - Gs * psi .* sN; F(:,:,:,1) F(:,:,:,1) w(i) * fcom * c(i,1); F(:,:,:,2) F(:,:,:,2) w(i) * fcom * c(i,2); F(:,:,:,3) F(:,:,:,3) w(i) * fcom * c(i,3); end end代码里circshift(psi, -c(i,:))取的是 xe_i 位置的伪势值因为circshift的第二个参数表示每个维度的移动量取负号才能拿到“下一个邻居”。fcom把流体-流体项和流体-固体项合并成一次累加方向分量分别写进三个通道。需要留意的是circshift隐含周期边界如果计算域四周不是固体或不是真正的周期边界力会串到对面。后面 4.3 节会给出一种掩码修正方案。提示指数伪势要求初始密度尽量接近 1。如果直接把初始密度设成几或者几十exp(-rho)的梯度会非常陡G 的稳定区间完全失效。这是 3D 版本最容易踩的坑。2.3 力怎么装进演化方程平衡态速度修正伪势力算出来后要进 LBM 的碰撞步。常见有两类接法一类是 Shan-Chen 早期论文里的平衡态速度修正另一类是 Guo 力模型。Guo 力模型把力作为离散源项加进演化方程对凝聚态伪势模型的控制更精细而平衡态速度修正实现简单教学代码里出现频率更高。用后者时先由密度场直接求宏观速度 u再把力折算进平衡态速度u_eq u τ F / ρ然后把 u_eq 代入平衡态分布函数。这样做的物理含义是外力先改变粒子的动量再参与碰撞松弛。碰撞步的向量化写法可以套 19 次循环没必要对每个网格点展开for i 1:19 cu c(i,1)*ueq(:,:,:,1) c(i,2)*ueq(:,:,:,2) c(i,3)*ueq(:,:,:,3); u2 ueq(:,:,:,1).^2 ueq(:,:,:,2).^2 ueq(:,:,:,3).^2; feq(:,:,:,i) w(i)*rho .* (1 3*cu 4.5*cu.^2 - 1.5*u2); f(:,:,:,i) f(:,:,:,i) - (f(:,:,:,i) - feq(:,:,:,i)) / tau; end注意这里声速平方取 1/3所以一阶项系数是 3二阶项系数是 4.5。τ 和运动粘度满足 ν (τ-0.5)/3格子单位因此 τ 越接近 0.5数值越容易震荡多孔介质里的常用区间是 0.7 到 1.0。G 决定两相密度差和表面张力但不能单独调得过大否则会把伪速度放大到和宏观驱替速度一个量级。下一章把这些代码串成完整的时间步。3. 在 MATLAB 中搭出 3D Shan-Chen 求解器骨架3.1 数据结构与初始条件三维 LBM 的分布函数存储常见格式是f(Nx,Ny,Nz,19)19 放最后一维。这样密度可以直接用sum(f,4)求19 个方向上的循环也不会破坏空间维度的内存布局。固体标记用逻辑数组solid(Nx,Ny,Nz)。固体节点上的分布函数依然需要占位否则计算宏观量时会出现 NaN处理方式是初始化为平衡态分布但后续碰撞和流步跳过这些节点。初始条件下面的写法比较稳Nx 32; Ny 32; Nz 32; rho0 ones(Nx,Ny,Nz); rng(0); rho0 rho0 0.001*randn(Nx,Ny,Nz); % 低幅扰动 f zeros(Nx,Ny,Nz,19); [c,w] d3q19(); for i 1:19 f(:,:,:,i) w(i)*rho0; end随机扰动的幅值要小因为它只负责打破对称性。Shan-Chen 模型在 G 足够大时会自己把密度拉成高密度相和低密度相初始不需要人工画出气泡。如果扰动幅值过大相当于给系统一个强力初始速度容易在第一轮时间步就发散。3.2 时间步碰撞、流步、反弹整个时间步用函数封装方便后面加输出和保存。流步用节点循环作为教学版本逻辑最清晰function [f, rho, u] lbm3d_step(f, solid, G, Gs, tau) [c,w] d3q19(); opp zeros(19,1); for i 1:19 [~,opp(i)] ismember(-c(i,:), c, rows); end rho sum(f,4); u zeros([size(rho),3]); for i 1:19 u(:,:,:,1) u(:,:,:,1) c(i,1)*f(:,:,:,i); u(:,:,:,2) u(:,:,:,2) c(i,2)*f(:,:,:,i); u(:,:,:,3) u(:,:,:,3) c(i,3)*f(:,:,:,i); end for d 1:3 u(:,:,:,d) u(:,:,:,d) ./ rho; end F shanchen_force(rho, solid, G, Gs); ueq u tau * F ./ max(rho, 1e-12); % 碰撞 for i 1:19 cu c(i,1)*ueq(:,:,:,1) c(i,2)*ueq(:,:,:,2) c(i,3)*ueq(:,:,:,3); u2 ueq(:,:,:,1).^2 ueq(:,:,:,2).^2 ueq(:,:,:,3).^2; feq w(i)*rho .* (1 3*cu 4.5*cu.^2 - 1.5*u2); f(:,:,:,i) f(:,:,:,i) - (f(:,:,:,i) - feq) / tau; end % 流步 半程反弹 [Nx,Ny,Nz] size(rho); fstream zeros(size(f)); for x 1:Nx for y 1:Ny for z 1:Nz if solid(x,y,z) continue; end for i 1:19 xi mod(x c(i,1) - 1, Nx) 1; yi mod(y c(i,2) - 1, Ny) 1; zi mod(z c(i,3) - 1, Nz) 1; if solid(xi,yi,zi) fstream(x,y,z,opp(i)) fstream(x,y,z,opp(i)) f(x,y,z,i); else fstream(xi,yi,zi,i) fstream(xi,yi,zi,i) f(x,y,z,i); end end end end end f fstream; rho sum(f,4); end这个版本把“目标节点是固体”的分布函数直接反弹回原格点等价于半程反弹能自动满足多孔介质内表面的无滑移条件。opp数组在函数开头用ismember匹配反向速度索引避免在多层循环里一次次搜索。ueq计算时用max(rho,1e-12)防止密度为零的网格除零。这套代码在 32 的三次方网格上跑几百步没有问题再大会很慢性能优化留到第 5 章。主循环只需要对 f、固体标记和两个控制参数做简单迭代steps 1500; for step 1:steps [f, rho, u] lbm3d_step(f, solid, 1.0, -0.05, 0.8); if mod(step, 100) 0 fprintf(step%d maxU%.3e mass%.6f\n, ... step, max(abs(u(:))), sum(rho(:))); end end其中 G 取 1.0 通常能形成明显两相分离Gs 取负的 -0.05 会削弱固体对非湿相的排斥具体符号取决于你希望哪一相挂在固体表面。masssum(rho(:))是每一轮都应该盯着的量它长时间明显漂移说明边界或力项有泄漏。3.3 常见运行错误怎么判断三维修正后的常见症状有两类。第一类是前几十步就出现 NaN绝大多数原因是 τ 太接近 0.5或者初始扰动太大少数情况是力项里出现了负密度。第二类是密度场不分离而是整体成噪声状通常是 G 太小界面张力无法克服伪速度扰动。可以把 G 从 0.6 开始每隔 0.1 往上试每次只跑 200 步看是否出现两个明显密度峰。高质量的两相状态密度直方图应该呈现双峰而不是一个拖尾单峰。4. 构造 3D 多孔介质并初始化两相驱替4.1 随机球堆叠生成孔隙骨架多孔介质数字化最常用的方法是在计算域内随机放置若干球体把球体内部标记为固体。这个方法能快速得到连通孔隙、死端孔和局部窄喉足够用来验证 Shan-Chen 模型的行为。function solid make_porous_random_spheres(N, r, nsphere) solid false(N,N,N); for k 1:nsphere cx randi([r1, N-r]); cy randi([r1, N-r]); cz randi([r1, N-r]); [X,Y,Z] ndgrid(1:N,1:N,1:N); sphereMask (X-cx).^2 (Y-cy).^2 (Z-cz).^2 r^2; solid solid | sphereMask; end porosity 1 - sum(solid(:)) / N^3; fprintf(porosity %.3f\n, porosity); end球体之间允许重叠重叠会让实际孔隙度比1-nsphere*(4/3)pi r^3/N^3高。所以要按目标孔隙度反推球数时需要跑两三组采样取平均。球体边界留了半径 r 的距离避免球被周期边界截断。ndgrid对 N48 和几十个球完全够用如果 N 超过 128内存会明显吃紧改成按球心局部填充更快。真正的岩心数据也可以用同样格式接进来。CT 扫描得到的三维二值掩码只要满足“solid 为 true 表示不可流动”就能直接替换这个函数。不同来源的孔隙度差异很大但 Shan-Chen 模型本身不依赖孔隙骨架是怎么生成的只依赖固体标记这也是 LBM 做多孔介质的一大便利。4.2 初始密度场和驱替相种子两相驱替需要给系统一个非平衡的初始条件否则两相在静态下只是均匀混合后相分离不会形成宏观定向流动。常见做法是在入口区域放置一个密度略高的种子相让伪势力在局部产生压力梯度之后入口边界持续补充该相。教学版本里最稳的不是直接设 2.0 这类大密度而是用扰动加区域种子rho ones(N,N,N) 0.001*randn(N,N,N); % 入口一个圆柱区域作为被驱替相种子 [X,Y,Z] ndgrid(1:N,1:N,1:N); seed (X 4) ((Y-N/2).^2 (Z-N/2).^2 6^2); rho(seed) 1.3;种子密度选择 1.3 而不是更高是为了避免界面压力突变带来的速射扰动。驱替方向沿着 x 轴入口截面上的高密度区域会逐渐向孔隙内推进。如果希望模拟润湿相驱替非湿相就把 Gs 调整为负值让密度较高的一相更愿意贴近固体表面反之取正值。4.3 力计算里的周期边界泄漏问题2.2 节留下的隐患现在要处理。circshift会让边界的邻居取到对面网格这对球体骨架内部没有影响但计算域外表面若没有固体力场会在外边界产生一条虚假的“界面”。最省事的解决方法是把计算域六个外表面全部标记为固体solid([1 end], :, :) true; solid(:, [1 end], :) true; solid(:, :, [1 end]) true;加上这一层后circshift在边界取到的邻居是固体伪势力只和流体-固体项有关不再和对面流体发生作用。配合流步里的反弹计算域外边界变成封闭壁面适合模拟静态接触角测试。如果要模拟流动还需要把入口和出口的固体标记去掉换成压力边界或者速度边界这部分通常用 Zou-He 边界实现教学实现可以先用周期边界跑通驱替后再换出入口边界。4.4 怎么看驱替结果三维数据直接可视化容易看不清内部孔隙。常见做法是取中间切片把固体和两相分开展示mid round(N/2); imagesc(squeeze(rho(:,mid,:))); axis equal; colormap(flipud(gray)); hold on; contour(squeeze(solid(:,mid,:)), [0.5 0.5], r);密度切片能看到两相界面红色轮廓对应固体骨架。另外还要统计平均密度和相饱和度相饱和度可以把密度阈值设成全局平均密度来分相。每 100 步计算一次检查驱替前缘是否在移动、饱和度是否单调变化。如果饱和度几百步都不变化先怀疑 G 太小导致驱替压力不够再检查入口种子是否被压力波冲散。5. Shan-Chen 参数稳定区间、3D 加速与结果验证5.1 调参顺序和参数表3D 多孔介质的 Shan-Chen 参数调起来比二维修长之后棘手因为界面的伪速度会被三维拓扑放大。我一般按固定顺序调参避免同时改两个变量导致问题定位困难。参数常用范围作用τ0.7 ~ 1.0粘度与稳定性越小越容易震荡G0.6 ~ 1.2两相密度差和表面张力过大产生高伪速度Gs-0.2 ~ 0.2固体润湿性负值亲液rho00.9 ~ 1.1初始密度基准配合指数伪势先固定 τ0.8、Gs0调 G 到两相分离再固定 G调 τ 让速度场平滑最后用 Gs 调接触角。每一轮只看一个指标比如最大宏观速度max(abs(u(:)))。驱替过程中这个值保持在 0.01 到 0.05 之间比较安全超过 0.1 就很容易在窄孔隙处产生局部震荡。5.2 让 MATLAB 版本跑得更快教学循环在 32³ 网格上没问题到 64³ 或 128³ 就扛不住了。常见做法是把流步的节点循环改写为预索引数组避免在时间步里反复计算目标坐标。先算好每个方向的目标线性索引再用f(srcIndex)一次性迁移。这个思路比circshift灵活也能兼容非周期边界ind zeros(Nx,Ny,Nz,19); for i 1:19 xi mod((1:Nx) c(i,1) - 1, Nx) 1; yi mod((1:Ny) c(i,2) - 1, Ny) 1; zi mod((1:Nz) c(i,3) - 1, Nz) 1; ind(:,:,:,i) sub2ind([Nx Ny Nz], ... repmat(xi,1,Ny,Nz), repmat(yi,Nx,1,Nz), repmat(zi,Nx,Ny,1)); end迁移时对每个方向做一次线性索引赋值反弹部分再对固体邻居做掩码修正。这样分布函数数组可以整体移动省掉三层循环。如果电脑有 NVIDIA GPUgpuArray可以把碰撞步和力项直接搬到 GPU 上circshift和四维数组运算在 GPU 上也支持。但流步的散射类操作在 GPU 上便宜有限迁移加反弹整体加速通常在 5 倍左右不要指望 MATLAB 的 GPU 版本能跑出 CUDA 手写内核的性能。5.3 结果验证的三个指标Shan-Chen 代码跑起来后第一步做静态液滴验证。在纯周期域放一个球形高密度区域让系统达到稳态测量液滴半径和内外密度差再换不同 G 多跑几次观察是否满足 Laplace 定律。接触角验证则需要一个平板固体表面把 Gs 从负到正扫描量出接触角变化是否单调。第二步做质量守恒检查。理想情况下总质量在每个时间步的变化应小于 0.1%。检查方法在 3.2 节已经写在主循环里。如果总质量单调下降说明边界处理丢了分布函数如果上升则多半是力项和碰撞步之间重复累加了质量。伪势力本身不改变密度总和问题只会出在流步和边界上。第三步检查驱替后的剩余相分布是否合理。把多孔介质骨架画出来看剩余液滴是不是都集中在死角和小喉道附近如果剩余相均匀分布在所有孔隙表面说明润湿性参数没有起到区分作用先回 5.1 重新调 Gs。本文还有配套的精品资源点击获取