简介Shan-Chen-in-3D.rar为一份基于Lattice Boltzmann MethodLBM的3D多孔介质流动模拟MATLAB实现核心代码采用经典的Shan-Chen模型面向流体力学数值模拟、孔隙介质渗流分析及LBM算法学习者。该模型通过势能函数模拟流体间的相互作用力可处理多相流、多组分流等非均匀流体问题在土壤、岩石等孔隙结构中的渗透率、流速分布和压力梯度计算方面具有典型应用对地下水资源管理、地质储油和环境修复等领域具有重要参考意义。压缩包内仅包含1个m格式源文件文件大小约2KB结构简明便于直接导入MATLAB运行与断点调试借助相关注释与输出结果可快速验证三维流动模拟的正确性。目前已有273人学习下载无论您是刚开始接触LBM的研究生还是从事流动仿真与多孔介质研究的工程师都能从中获得实用的算法示例与建模思路。通过逐行阅读并拆解这段代码读者可深入理解Shan-Chen模型在三维空间中的初始化、相互作用力计算、碰撞迁移等核心步骤积累多孔介质边界条件与参数设置的实战经验并为后续使用OpenLB等开源LBM平台开展更复杂模拟提供直接对照与扩展基础。1. 三维Shan-Chen多孔介质LBM直接从文件名到驱动非混相流的MATLAB骨架“Shan-Chen-in-3D.rar_3D porous_LBM Shan-Chen_LBM shan_matlab shan”这个压缩包放到实际需求里通常不是“跑通一个演示动画”而是在三维多孔结构里把两相驱替算到能出毛细压力曲线、相对渗透率曲线、残余油饱和度的程度。LBM里能扛这个任务的多数人第一选择就是Shan-Chen伪势模型它不用显式追踪相界面每个格点算一个伪势力用状态方程参数控制密度比和表面张力实现代码短、参数直观MATLAB的矩阵索引和切片可视化又天然适合逐层观察界面推进。这套骨架能解决的具体问题很明确给一张二值化的孔隙骨架、两组流体属性输出谁先占据哪条孔道、何时发生卡断、不同毛管数下的相分布。对新手它能带你进入三维两相LBM对老手它能作为换成MRT或GPU时参照的基线版本。2. 三维Shan-Chen模型的前置伪势力、D3Q19速度集与外力接入方式2.1 为什么三维不能照搬二维的D2Q9二维LBM常用D2Q9九个速度方向在三维里不够用z方向的动量传递无法表达固体格点的拓扑连通也会失真。三维Shan-Chen实现必须至少用D3Q19它包含1个静止方向、6个面心和12个对角方向能覆盖三维空间中所有轴对齐和面对角传播。少数代码用D3Q15或D3Q27前者各向同性略差后者方向多但计算量上涨明显。对多孔介质驱替这种“统计平均量比瞬时细节更重要”的场景D3Q19是性价比最稳的配置。对应的权重系数表是每个方向求和时的分母直接决定离散拉普拉斯和各向同性误差。实际编码时我把19个方向按“静止、轴向、对角”三类分组存成稀疏表% D3Q19 速度分量每行一个方向c [cx cy cz] c zeros(19,3); c(1,:) [0 0 0]; % 轴向 6 个方向 c(2:7,:) [1 0 0; -1 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 -1]; % 对角 12 个方向每行包含两个非零分量 c(8:19,:) [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 zeros(19,1); w(1) 1/3; w(2:7) 1/18; w(8:19) 1/36;这里关键是别把权重记成“1/18和1/36都是对的但顺序要和c矩阵行号对齐”。碰撞和迁移循环里权重是按下标直接取用的顺序错位会让力的方向各向异性表现为静止液滴压强场呈方形而非圆形。多数新手调不出球形液滴第一步错的就在这张表。2.2 Shan-Chen伪势力在三维网格上的梯度写法SC模型的核心是伪势函数和局域相互作用力。常用伪势形式是ψ(ρ) 1 - exp(-ρ)它的好处是当ρ增大时ψ渐近到1数值上稳定不会因为密度峰值过高导致力发散。三维多孔介质里密度分布极不均匀强疏水壁面附近流体密度会剧烈变化用这个截断形式能避免局部“闪蒸”伪影。力项的标准写法是F(x,t) -G ψ(x) Σq wq ψ(x cq) cq其中求和遍历邻近19个方向。这段代码用MATLAB的circshift实现周期邻居索引是三维实现里最干净的写法function [Fx,Fy,Fz] sc_force3d(rho, G, psi_weights, c) % rho: 三维密度场 % G: 相互作用强度负值代表分子间引力 psi 1 - exp(-rho); [Fx,Fy,Fz] deal(zeros(size(rho))); for q 1:19 if c(q,1)0 c(q,2)0 c(q,3)0 continue; end % circshift 负数表示向该方向取邻居 psi_shift circshift(psi, -c(q,:)); Fx Fx psi_weights(q) * psi_shift * c(q,1); Fy Fy psi_weights(q) * psi_shift * c(q,2); Fz Fz psi_weights(q) * psi_shift * c(q,3); end Fx -G .* psi .* Fx; Fy -G .* psi .* Fy; Fz -G .* psi .* Fz; % 归一化防止在锐利界面处力过大导致发散 maxF max(sqrt(Fx.^2 Fy.^2 Fz.^2), [], all); if maxF 0.02 scale 0.02 / maxF; Fx Fx * scale; Fy Fy * scale; Fz Fz * scale; end end这个力项的计算量不小每个格点要做18次邻居读取和累加。MATLAB里直接全网格算三维网格到128³就明显变慢常见做法是预先建一个流体格点索引表只对流体格点算力固体格点跳过。上面代码的归一化阈值不是LBM理论里的固定参数但实际调试中强烈建议留这条线最多牺牲一点表面张力量值换来计算稳定性。2.3 外力入演化方程Guo格式与速度偏移计算出的F不是直接加进分布函数。LBM演化方程里外力项和宏观速度之间存在交叉项跳过它会导致渗透率系统偏差。常用做法是Guo格式碰撞前先用含外力的速度重定义平衡态速度。% 在碰撞-迁移主循环里 omega 1/tau; u (ux Fx./(2*rho)); v (uy Fy./(2*rho)); we (uz Fz./(2*rho)); % 平衡态分布函数D3Q19公式 for q 1:19 cu c(q,1)*u c(q,2)*v c(q,3)*we; feq(:,:,:,q) w(q) * rho .* (1 3*cu 4.5*cu.^2 - 1.5*(u.^2v.^2we.^2)); end % 碰撞含外力 f f omega * (feq - f) fc_guo;外力项加在碰撞步的末尾而不是直接改分布函数值。这样处理能保证连续性方程在离散层面仍然守恒驱替过程中的质量损失会小一个量级。特别注意这里的速度要减掉“半力修正”否则宏观速度会偏大渗透率模拟值漂移约0.5个tau系数。3. 在MATLAB里生成真实可用的三维多孔骨架体素标记、连通域剪枝与周期性边界3.1 用随机球体堆叠生成带可调孔隙度的骨架多孔介质LBM的第一件事不是写演化代码而是先把固体骨架制造出来。最常用的生成方法是在三维网格里随机放球体球心坐标服从均匀分布半径可以固定也可以加随机扰动。这个模型接近真实砂粒堆积而且生成速度快不需要CT扫描数据也能支撑算法验证。Nx 96; Ny 96; Nz 96; solid false(Nx, Ny, Nz); rng(2026); % 固定随机种子结果可复现 num_seeds 400; % 种子数越多孔隙度越低 R_mean 6; % 球体平均半径单位格点 [Xg, Yg, Zg] ndgrid(1:Nx, 1:Ny, 1:Nz); for k 1:num_seeds cx randi(Nx); cy randi(Ny); cz randi(Nz); R R_mean (rand-0.5)*2; % 半径扰动 dist2 (Xg-cx).^2 (Yg-cy).^2 (Zg-cz).^2; solid(dist2 R^2) true; end region ~solid; porosity sum(region(:)) / numel(region); fprintf(当前孔隙度: %.3f\n, porosity);数种子数是一次性标定的通常用孔隙度反推目标孔隙度0.3、球半径6在96³网格里约需400到600个种子。每次生成后打印孔隙度偏了调整种子数再跑不要指望一次到位。这个代码段的坑在于ndgrid生成的是Mx×My×Mz的网格矩阵距离计算会占很大内存96³已经接近单机普通内存的舒适区更大网格建议按层生成。3.2 用连通域分析剪掉死端孔随机球体堆叠会附带许多孤立的孔隙团簇这些不连通区域在真实岩心里大概率是盲端或完全封闭的。驱替模拟若不过滤这些区域入口压力根本传不进去流体只会聚集在连通主干外的小洞里结果把毛管压力算成负数。参数调试时最迷惑的现象——残余饱和度异常高——往往不是Shan-Chen参数问题而是连通域没处理。% 连通域分析6-邻域连通即共享一个面视为连通 cc bwconncomp(region, 6); num_clusters cc.NumObjects; pore_label labelmatrix(cc); % 取域中心附近的一个孔隙格点作为主连通域的标号 anchor pore_label(Nx/2, Ny/2, Nz/2); if anchor 0 error(中心点被固体占据换一个anchor点); end main_region (pore_label anchor); fprintf(主连通域占全部孔隙的比例: %.3f\n, ... sum(main_region(:)) / sum(region(:)));如果主连通域占比小于85%说明生成的骨架质量不好应增大球半径或减少种子数。低于70%时模拟效率极低因为你在一堆无连接的空间里白白计算。过滤后得到的main_region才是后续LBM计算的孔隙掩膜。这个步骤看似多余却是多孔流模拟结果具备统计意义的前提。3.3 入口、出口和周期性边界的处理规范矿山级多孔介质模拟最头疼的是边界怎么设。全周期边界在数学上最干净但物理上你没法写“从左边注入、右边采出”。最常见的工程做法是实际骨架只放在中间区域入口面和出口面各留几层全孔隙缓冲区然后沿流动方向作密度驱动。% 沿z方向驱替预留入口出口缓冲区 buffer_z 6; solid(:, :, 1:buffer_z) 0; % 入口缓冲区全部置为孔隙 solid(:, :, end-buffer_z1:end) 0; % 出口缓冲区 % 固体标记0流体1固体 region ~solid main_region_valid; region(:, :, 1:buffer_z) 1; % 入口缓冲区强制为孔隙 region(:, :, end-buffer_z1:end) 1;这里的缓冲区同时承担“进出口压力条件”的载体。简单做法可以像文献里那样用恒定密度注入在入口缓冲区每一层固定设置驱替流体密度ρ_in出口固定为另一相密度ρ_out这样形成一个稳定压力梯度。缓冲区厚度建议最少6格太薄会产生边界数值反射界面到达边界时产生多余的振荡。3D模拟中这个边界误差会沿着x和y方向扩散比二维严重。4. 三维Shan-Chen的MATLAB实现分布函数数组、碰撞迁移循环与固体格点处理4.1 分布函数的存储与维度排列分布函数的存储方式直接决定性能。常见方案是把f定义成四维数组 f(Nx,Ny,Nz,19)第四个维度存19个方向。这个布局下碰撞步的向量化比较容易但迁移步要循环19次把每个方向的分量按速度向量移位。% f: 第4维放19个方向内存约 Nx*Ny*Nz*19*8 字节 % 96^3 网格double类型约 126MB用single可减半 f zeros(Nx, Ny, Nz, 19, single); rho zeros(Nx, Ny, Nz, single);对130³以上网格建议直接single。LBM的浮点误差主要是截断误差single精度对多孔介质这种统计平均量足够节省下来的内存能扩展一到两倍网格规模对分辨孔喉尤为重要。初始化时所有孔隙格点赋予平衡态分布固体格点f全部置零并在每个时间步维持不要参与碰撞。4.2 碰撞-迁移循环的主函数骨架主循环里三件事算宏观量、执行碰撞含外力修正、执行迁移。三维迁移每个方向都是移一层数组用circshift一行完成for it 1:max_steps % 1. 宏观量 rho sum(f, 4); % 2. 速度先算动量再除密度 ux zeros(Nx,Ny,Nz,single); uy ux; uz ux; for q 1:19 ux ux f(:,:,:,q) * c(q,1); uy uy f(:,:,:,q) * c(q,2); uz uz f(:,:,:,q) * c(q,3); end % 3. Shan-Chen力只对孔隙区计算 [fx, fy, fz] sc_force3d(rho, G, w, c); % 4. 将外力引入速度 ux (ux 0.5*fx) ./ rho; uy (uy 0.5*fy) ./ rho; uz (uz 0.5*fz) ./ rho; % 5. 碰撞 % 遍历方向构造平衡态 for q 1:19 cu c(q,1)*ux c(q,2)*uy c(q,3)*uz; u2 ux.^2 uy.^2 uz.^2; feq w(q) * rho .* (1 3*cu 4.5*cu.^2 - 1.5*u2); f(:,:,:,q) f(:,:,:,q) omega * (feq - f(:,:,:,q)); end % 6. 迁移每个方向把对应分量移位 f_shifted zeros(size(f), single); for q 1:19 f_shifted(:,:,:,q) circshift(f(:,:,:,q), -c(q,:)); end f f_shifted; % 7. 固体格点反弹边界处理 % 固体格点反射把进入固体的分布按反方向弹回相邻流体格点 [fx_solid, fy_solid, fz_solid] ... % 简化示意 % 真实实现对固体点逐方向做 f(opp) f(neighbor) end这个骨架是所有后续工作的地基。第4步的速度修正来自Guo格式如果省掉“0.5*f”宏观速度和实际动量之间会产生系统性偏移表现为渗透率偏高约3%到5%。第7步的边界条件用“半步反弹”最常用当某个方向的分布函数迁移到固体格点时把它弹回原流体格点的反向方向。在三维多孔里固体点数远多于流体邻接数遍历固体mask每点做一次弹回即可。4.3 空孔隙区域与多相初始池布置启动两相驱替前要在连通孔隙里放置初始流体。简化做法是把连通域按z方向切两半下半填充湿相例如水上半填充非湿相例如油中间留几个格点的过渡带让界面自然松弛。这个初始化会立刻触发高密度梯度所以头两三百步要用松弛时间τ1.0以上的大粘度作初值逐步降到目标值。这个“预热”阶段不算物理结果只用来让界面光滑过渡。5. 核心参数标定G、ρ0、τ、Gads这4个必调量怎么和物理对应5.1 状态方程强度G与密度比的关系Shan-Chen模型里G不是随意设的。G的绝对值越大两相密度差越大但超过某个数值域就会触发数值失稳出现负密度或周期性振荡。经验做法是固定伪势形式ψ 1-exp(-ρ)、平均密度ρ0 1.0时G的绝对值设在0.1到0.2之间。标定密度比的方法是先在均匀网格里放一个小液滴平衡后统计液滴内部密度ρ_l和外部密度ρ_v打印出来看是否满足目标比。% 在无固体域做液滴标定 G -0.15; rho0 1.0; rho repmat(rho0, [Nx, Ny, Nz]); % 在中心放一个半径R15的球初始密度为外界1.5倍 [Xg,Yg,Zg] ndgrid(1:Nx,1:Ny,1:Nz); inner (Xg-Nx/2).^2 (Yg-Ny/2).^2 (Zg-Nz/2).^2 15^2; rho(inner) 1.5; % 运行 SC 演化 2000 步 % 统计流量 rl mean(rho(inner rho0.4)); rv mean(rho(~inner rho0.4)); fprintf(密度比: %.2f\n, rl/rv);标定时要确认平衡后的密度不是初始给的1.5和0.8而是SC状态方程自己收敛出的值。密度比一旦低于目标档位加大G即可不用动其他参数。但G绝对值超过0.2后力尖峰会变得不可收敛此时应该换用更陡的伪势形如ψ ρ * exp(-ρ/ρ0)的变体工程上也比较常见。5.2 松弛时间τ与粘度、渗透率的换算τ决定流体运动粘度ν cs²(τ - 0.5)其中cs²1/3。多孔介质模拟对τ非常敏感τ太大边界层变厚微小孔隙内的流动被过度阻尼τ太接近0.5数值振荡增大。工程区间在0.7到1.5之间我自己一般固定在τ1.0起步验证通过后再改。用已知渗透率的解析阵列比如一排直圆柱组成的规则阵列标定模型渗透率时公式为K ν · Q · L / (Δt · ΔP)其中Q是体积流量、L是渗流长度。如果标定结果比解析值偏差超过5%先检查入口/出口缓冲区和力项速度修正这两个因素占误差来源大头。5.3 壁面吸附力Gads与接触角的关系接触角由固壁对两相的相对亲和性决定。Shan-Chen里常见的做法是在壁面格点施加额外吸附力Gads壁面处流体所受额外力与Gads符号相关Gads为正时壁面吸引流体中密度较高的部分接触角变小表现为亲液。标定接触角的标准手段是在一个平板通道内放一个静态液滴平衡后量三相接触点上的切向夹角。% 平板接触角标定z方向上下壁面 Gads 0.06; % 正亲湿相负憎湿相 % 在作用力函数里固体mask为1处额外加 for q 1:19 % 如果邻居是固体则把壁面密度当作常数rho_wall rho_wall 1.0; psi_w 1 - exp(-rho_wall); Fx Fx Gads * psi .* psi_w .* w(q) * c(q,1); end投影接触角与Gads的关系不是线性的通常要做三到四个不同Gads的标定实验拟合Gads-cosθ关系曲线。这样你在指定目标接触角后可以直接反查该用哪个Gads而不是盲猜之后反复跑全尺寸模拟。固体壁面的“刚度”也要注意如果壁面伪势太小界面临近壁面时会提前失稳出现不该有的附着层。5.4 时间步长与收敛判定SC模型的时间步长在LBM里是隐含的每个时间步对应物理时间Δt多用无量纲毛管数Ca来控制。多孔介质驱替里Ca μu/σ设目标Ca1e-4到1e-3反算驱替速度或者流量。运行时判断收敛不是看某个参数而是观察出口流量随时间是否进入平台期。这里给出常见的参数参考表表格里的量纲均为格子单位参数符号常用范围典型初值相互作用强度G-0.20 ~ -0.10-0.15平均密度ρ00.8 ~ 2.01.0松弛时间τ0.7 ~ 1.51.0壁面吸附力Gads-0.20 ~ 0.200.06伪势截断系数无1.0 ~ 2.01.0这些参数不是独立调节的。动了G之后密度比变化界面厚度也跟着变接触角标定会失准所以业界经验是先固定G调到目标密度比再固定τ调到目标粘度最后调Gads配接触角。中途换G前面标定全部重来。6. 验证与排错用Laplace标定法确认三维Shan-Chen实现可靠6.1 两球法测表面张力的最小代码实现SC之后第一件事不是跑多孔而是先在空旷网格里验证界面张力是否满足Laplace定律。做法是在网格中心放一个球形液滴等系统平衡后测液滴内外压差然后由ΔP 2σ/R反算σ。如果σ随R变化明显偏离常数说明力项离散有问题。% 静态液滴验证 R_list [12, 16, 20]; sigma_list zeros(3,1); for i 1:3 rho ones(Nx,Ny,Nz,single); R R_list(i); [Xg,Yg,Zg] ndgrid(1:Nx,1:Ny,1:Nz); inner (Xg-Nx/2).^2 (Yg-Ny/2).^2 (Zg-Nz/2).^2 R^2; rho(inner) 1.5; % 液滴初始密度 run_evolution(2000); % 取两相区域的体积平均压力 p_l mean(rho(inner) * rho0 / 3); p_v mean(rho(~inner) * rho0 / 3); sigma_list(i) (p_l - p_v) * R / 2; end if std(sigma_list) / mean(sigma_list) 0.05 fprintf(Laplace 验证通过σ≈%.4f\n, mean(sigma_list)); else warning(σ波动过大检查速度集和力项权重); end这个测试跑不通后面多孔模拟的毛管力全不可信。常见失败点是压力统计取了包含界面过渡层的格点统计区域要离界面至少两个格子否则压差偏大。6.2 三个高频排错点第一是质量守恒异常。每500步输出一次总质量sum(rho)若随时间单调增长或衰减超过0.5%要么是边界泄漏要么是Guo外力项在固体附近处理错误。第二是界面“碎成雾状”的抖动典型原因是G过大或τ小于0.7拉回参数区间即可如果是G值没问题但界面仍然抖动就检查伪势ψ是否在某些格点出现负值这种情况多半是密度初始扰动过猛。第三是驱替完全没有进展先看入口缓冲区是否被固体堵住再看连通域是否有中心孤岛导致压差传不进去。6.3 三维界面观察的一个实用技巧最后给一个我常用的MATLAB可视化检查法——不等时间步就调参数而是把驱替界面保存成isosurface数据若干步后再画。把f和rho的view只保留关键等值面会让内存占用小很多这一步用MATLAB自带的isosurface和patch画组相界面位置能直接看出前沿是否干净、是否有孤立液滴被拖后。% 用等值面提取两相界面 lv 0.5; % 密度等值面阈值两相密度中间值附近 [fv, fvc] isosurface(Xg, Yg, Zg, rho, lv); p patch(Faces, fv, Vertices, fvc); p.FaceColor interp; p.EdgeColor none; view(3); camlight; xlabel(x); ylabel(y); zlabel(z);这里核心技巧是不要直接画三维体数据体素太大MATLAB画图会卡到没法交互只用isosurface提取界面轮廓。界面形态从光滑曲面变成锯齿状碎片说明当前参数接近数值失稳边界及时降G或升高ρ0能救活很多本要重跑的长任务。参数没问题时这个界面应该会把孔隙通道逐步切开而不是在孔喉处留一条细长的灰色拖尾。本文还有配套的精品资源点击获取 三维Shan-Chen多孔介质LBM的MATLAB实现与参数标定指南
简介Shan-Chen-in-3D.rar为一份基于Lattice Boltzmann MethodLBM的3D多孔介质流动模拟MATLAB实现核心代码采用经典的Shan-Chen模型面向流体力学数值模拟、孔隙介质渗流分析及LBM算法学习者。该模型通过势能函数模拟流体间的相互作用力可处理多相流、多组分流等非均匀流体问题在土壤、岩石等孔隙结构中的渗透率、流速分布和压力梯度计算方面具有典型应用对地下水资源管理、地质储油和环境修复等领域具有重要参考意义。压缩包内仅包含1个m格式源文件文件大小约2KB结构简明便于直接导入MATLAB运行与断点调试借助相关注释与输出结果可快速验证三维流动模拟的正确性。目前已有273人学习下载无论您是刚开始接触LBM的研究生还是从事流动仿真与多孔介质研究的工程师都能从中获得实用的算法示例与建模思路。通过逐行阅读并拆解这段代码读者可深入理解Shan-Chen模型在三维空间中的初始化、相互作用力计算、碰撞迁移等核心步骤积累多孔介质边界条件与参数设置的实战经验并为后续使用OpenLB等开源LBM平台开展更复杂模拟提供直接对照与扩展基础。1. 三维Shan-Chen多孔介质LBM直接从文件名到驱动非混相流的MATLAB骨架“Shan-Chen-in-3D.rar_3D porous_LBM Shan-Chen_LBM shan_matlab shan”这个压缩包放到实际需求里通常不是“跑通一个演示动画”而是在三维多孔结构里把两相驱替算到能出毛细压力曲线、相对渗透率曲线、残余油饱和度的程度。LBM里能扛这个任务的多数人第一选择就是Shan-Chen伪势模型它不用显式追踪相界面每个格点算一个伪势力用状态方程参数控制密度比和表面张力实现代码短、参数直观MATLAB的矩阵索引和切片可视化又天然适合逐层观察界面推进。这套骨架能解决的具体问题很明确给一张二值化的孔隙骨架、两组流体属性输出谁先占据哪条孔道、何时发生卡断、不同毛管数下的相分布。对新手它能带你进入三维两相LBM对老手它能作为换成MRT或GPU时参照的基线版本。2. 三维Shan-Chen模型的前置伪势力、D3Q19速度集与外力接入方式2.1 为什么三维不能照搬二维的D2Q9二维LBM常用D2Q9九个速度方向在三维里不够用z方向的动量传递无法表达固体格点的拓扑连通也会失真。三维Shan-Chen实现必须至少用D3Q19它包含1个静止方向、6个面心和12个对角方向能覆盖三维空间中所有轴对齐和面对角传播。少数代码用D3Q15或D3Q27前者各向同性略差后者方向多但计算量上涨明显。对多孔介质驱替这种“统计平均量比瞬时细节更重要”的场景D3Q19是性价比最稳的配置。对应的权重系数表是每个方向求和时的分母直接决定离散拉普拉斯和各向同性误差。实际编码时我把19个方向按“静止、轴向、对角”三类分组存成稀疏表% D3Q19 速度分量每行一个方向c [cx cy cz] c zeros(19,3); c(1,:) [0 0 0]; % 轴向 6 个方向 c(2:7,:) [1 0 0; -1 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 -1]; % 对角 12 个方向每行包含两个非零分量 c(8:19,:) [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 zeros(19,1); w(1) 1/3; w(2:7) 1/18; w(8:19) 1/36;这里关键是别把权重记成“1/18和1/36都是对的但顺序要和c矩阵行号对齐”。碰撞和迁移循环里权重是按下标直接取用的顺序错位会让力的方向各向异性表现为静止液滴压强场呈方形而非圆形。多数新手调不出球形液滴第一步错的就在这张表。2.2 Shan-Chen伪势力在三维网格上的梯度写法SC模型的核心是伪势函数和局域相互作用力。常用伪势形式是ψ(ρ) 1 - exp(-ρ)它的好处是当ρ增大时ψ渐近到1数值上稳定不会因为密度峰值过高导致力发散。三维多孔介质里密度分布极不均匀强疏水壁面附近流体密度会剧烈变化用这个截断形式能避免局部“闪蒸”伪影。力项的标准写法是F(x,t) -G ψ(x) Σq wq ψ(x cq) cq其中求和遍历邻近19个方向。这段代码用MATLAB的circshift实现周期邻居索引是三维实现里最干净的写法function [Fx,Fy,Fz] sc_force3d(rho, G, psi_weights, c) % rho: 三维密度场 % G: 相互作用强度负值代表分子间引力 psi 1 - exp(-rho); [Fx,Fy,Fz] deal(zeros(size(rho))); for q 1:19 if c(q,1)0 c(q,2)0 c(q,3)0 continue; end % circshift 负数表示向该方向取邻居 psi_shift circshift(psi, -c(q,:)); Fx Fx psi_weights(q) * psi_shift * c(q,1); Fy Fy psi_weights(q) * psi_shift * c(q,2); Fz Fz psi_weights(q) * psi_shift * c(q,3); end Fx -G .* psi .* Fx; Fy -G .* psi .* Fy; Fz -G .* psi .* Fz; % 归一化防止在锐利界面处力过大导致发散 maxF max(sqrt(Fx.^2 Fy.^2 Fz.^2), [], all); if maxF 0.02 scale 0.02 / maxF; Fx Fx * scale; Fy Fy * scale; Fz Fz * scale; end end这个力项的计算量不小每个格点要做18次邻居读取和累加。MATLAB里直接全网格算三维网格到128³就明显变慢常见做法是预先建一个流体格点索引表只对流体格点算力固体格点跳过。上面代码的归一化阈值不是LBM理论里的固定参数但实际调试中强烈建议留这条线最多牺牲一点表面张力量值换来计算稳定性。2.3 外力入演化方程Guo格式与速度偏移计算出的F不是直接加进分布函数。LBM演化方程里外力项和宏观速度之间存在交叉项跳过它会导致渗透率系统偏差。常用做法是Guo格式碰撞前先用含外力的速度重定义平衡态速度。% 在碰撞-迁移主循环里 omega 1/tau; u (ux Fx./(2*rho)); v (uy Fy./(2*rho)); we (uz Fz./(2*rho)); % 平衡态分布函数D3Q19公式 for q 1:19 cu c(q,1)*u c(q,2)*v c(q,3)*we; feq(:,:,:,q) w(q) * rho .* (1 3*cu 4.5*cu.^2 - 1.5*(u.^2v.^2we.^2)); end % 碰撞含外力 f f omega * (feq - f) fc_guo;外力项加在碰撞步的末尾而不是直接改分布函数值。这样处理能保证连续性方程在离散层面仍然守恒驱替过程中的质量损失会小一个量级。特别注意这里的速度要减掉“半力修正”否则宏观速度会偏大渗透率模拟值漂移约0.5个tau系数。3. 在MATLAB里生成真实可用的三维多孔骨架体素标记、连通域剪枝与周期性边界3.1 用随机球体堆叠生成带可调孔隙度的骨架多孔介质LBM的第一件事不是写演化代码而是先把固体骨架制造出来。最常用的生成方法是在三维网格里随机放球体球心坐标服从均匀分布半径可以固定也可以加随机扰动。这个模型接近真实砂粒堆积而且生成速度快不需要CT扫描数据也能支撑算法验证。Nx 96; Ny 96; Nz 96; solid false(Nx, Ny, Nz); rng(2026); % 固定随机种子结果可复现 num_seeds 400; % 种子数越多孔隙度越低 R_mean 6; % 球体平均半径单位格点 [Xg, Yg, Zg] ndgrid(1:Nx, 1:Ny, 1:Nz); for k 1:num_seeds cx randi(Nx); cy randi(Ny); cz randi(Nz); R R_mean (rand-0.5)*2; % 半径扰动 dist2 (Xg-cx).^2 (Yg-cy).^2 (Zg-cz).^2; solid(dist2 R^2) true; end region ~solid; porosity sum(region(:)) / numel(region); fprintf(当前孔隙度: %.3f\n, porosity);数种子数是一次性标定的通常用孔隙度反推目标孔隙度0.3、球半径6在96³网格里约需400到600个种子。每次生成后打印孔隙度偏了调整种子数再跑不要指望一次到位。这个代码段的坑在于ndgrid生成的是Mx×My×Mz的网格矩阵距离计算会占很大内存96³已经接近单机普通内存的舒适区更大网格建议按层生成。3.2 用连通域分析剪掉死端孔随机球体堆叠会附带许多孤立的孔隙团簇这些不连通区域在真实岩心里大概率是盲端或完全封闭的。驱替模拟若不过滤这些区域入口压力根本传不进去流体只会聚集在连通主干外的小洞里结果把毛管压力算成负数。参数调试时最迷惑的现象——残余饱和度异常高——往往不是Shan-Chen参数问题而是连通域没处理。% 连通域分析6-邻域连通即共享一个面视为连通 cc bwconncomp(region, 6); num_clusters cc.NumObjects; pore_label labelmatrix(cc); % 取域中心附近的一个孔隙格点作为主连通域的标号 anchor pore_label(Nx/2, Ny/2, Nz/2); if anchor 0 error(中心点被固体占据换一个anchor点); end main_region (pore_label anchor); fprintf(主连通域占全部孔隙的比例: %.3f\n, ... sum(main_region(:)) / sum(region(:)));如果主连通域占比小于85%说明生成的骨架质量不好应增大球半径或减少种子数。低于70%时模拟效率极低因为你在一堆无连接的空间里白白计算。过滤后得到的main_region才是后续LBM计算的孔隙掩膜。这个步骤看似多余却是多孔流模拟结果具备统计意义的前提。3.3 入口、出口和周期性边界的处理规范矿山级多孔介质模拟最头疼的是边界怎么设。全周期边界在数学上最干净但物理上你没法写“从左边注入、右边采出”。最常见的工程做法是实际骨架只放在中间区域入口面和出口面各留几层全孔隙缓冲区然后沿流动方向作密度驱动。% 沿z方向驱替预留入口出口缓冲区 buffer_z 6; solid(:, :, 1:buffer_z) 0; % 入口缓冲区全部置为孔隙 solid(:, :, end-buffer_z1:end) 0; % 出口缓冲区 % 固体标记0流体1固体 region ~solid main_region_valid; region(:, :, 1:buffer_z) 1; % 入口缓冲区强制为孔隙 region(:, :, end-buffer_z1:end) 1;这里的缓冲区同时承担“进出口压力条件”的载体。简单做法可以像文献里那样用恒定密度注入在入口缓冲区每一层固定设置驱替流体密度ρ_in出口固定为另一相密度ρ_out这样形成一个稳定压力梯度。缓冲区厚度建议最少6格太薄会产生边界数值反射界面到达边界时产生多余的振荡。3D模拟中这个边界误差会沿着x和y方向扩散比二维严重。4. 三维Shan-Chen的MATLAB实现分布函数数组、碰撞迁移循环与固体格点处理4.1 分布函数的存储与维度排列分布函数的存储方式直接决定性能。常见方案是把f定义成四维数组 f(Nx,Ny,Nz,19)第四个维度存19个方向。这个布局下碰撞步的向量化比较容易但迁移步要循环19次把每个方向的分量按速度向量移位。% f: 第4维放19个方向内存约 Nx*Ny*Nz*19*8 字节 % 96^3 网格double类型约 126MB用single可减半 f zeros(Nx, Ny, Nz, 19, single); rho zeros(Nx, Ny, Nz, single);对130³以上网格建议直接single。LBM的浮点误差主要是截断误差single精度对多孔介质这种统计平均量足够节省下来的内存能扩展一到两倍网格规模对分辨孔喉尤为重要。初始化时所有孔隙格点赋予平衡态分布固体格点f全部置零并在每个时间步维持不要参与碰撞。4.2 碰撞-迁移循环的主函数骨架主循环里三件事算宏观量、执行碰撞含外力修正、执行迁移。三维迁移每个方向都是移一层数组用circshift一行完成for it 1:max_steps % 1. 宏观量 rho sum(f, 4); % 2. 速度先算动量再除密度 ux zeros(Nx,Ny,Nz,single); uy ux; uz ux; for q 1:19 ux ux f(:,:,:,q) * c(q,1); uy uy f(:,:,:,q) * c(q,2); uz uz f(:,:,:,q) * c(q,3); end % 3. Shan-Chen力只对孔隙区计算 [fx, fy, fz] sc_force3d(rho, G, w, c); % 4. 将外力引入速度 ux (ux 0.5*fx) ./ rho; uy (uy 0.5*fy) ./ rho; uz (uz 0.5*fz) ./ rho; % 5. 碰撞 % 遍历方向构造平衡态 for q 1:19 cu c(q,1)*ux c(q,2)*uy c(q,3)*uz; u2 ux.^2 uy.^2 uz.^2; feq w(q) * rho .* (1 3*cu 4.5*cu.^2 - 1.5*u2); f(:,:,:,q) f(:,:,:,q) omega * (feq - f(:,:,:,q)); end % 6. 迁移每个方向把对应分量移位 f_shifted zeros(size(f), single); for q 1:19 f_shifted(:,:,:,q) circshift(f(:,:,:,q), -c(q,:)); end f f_shifted; % 7. 固体格点反弹边界处理 % 固体格点反射把进入固体的分布按反方向弹回相邻流体格点 [fx_solid, fy_solid, fz_solid] ... % 简化示意 % 真实实现对固体点逐方向做 f(opp) f(neighbor) end这个骨架是所有后续工作的地基。第4步的速度修正来自Guo格式如果省掉“0.5*f”宏观速度和实际动量之间会产生系统性偏移表现为渗透率偏高约3%到5%。第7步的边界条件用“半步反弹”最常用当某个方向的分布函数迁移到固体格点时把它弹回原流体格点的反向方向。在三维多孔里固体点数远多于流体邻接数遍历固体mask每点做一次弹回即可。4.3 空孔隙区域与多相初始池布置启动两相驱替前要在连通孔隙里放置初始流体。简化做法是把连通域按z方向切两半下半填充湿相例如水上半填充非湿相例如油中间留几个格点的过渡带让界面自然松弛。这个初始化会立刻触发高密度梯度所以头两三百步要用松弛时间τ1.0以上的大粘度作初值逐步降到目标值。这个“预热”阶段不算物理结果只用来让界面光滑过渡。5. 核心参数标定G、ρ0、τ、Gads这4个必调量怎么和物理对应5.1 状态方程强度G与密度比的关系Shan-Chen模型里G不是随意设的。G的绝对值越大两相密度差越大但超过某个数值域就会触发数值失稳出现负密度或周期性振荡。经验做法是固定伪势形式ψ 1-exp(-ρ)、平均密度ρ0 1.0时G的绝对值设在0.1到0.2之间。标定密度比的方法是先在均匀网格里放一个小液滴平衡后统计液滴内部密度ρ_l和外部密度ρ_v打印出来看是否满足目标比。% 在无固体域做液滴标定 G -0.15; rho0 1.0; rho repmat(rho0, [Nx, Ny, Nz]); % 在中心放一个半径R15的球初始密度为外界1.5倍 [Xg,Yg,Zg] ndgrid(1:Nx,1:Ny,1:Nz); inner (Xg-Nx/2).^2 (Yg-Ny/2).^2 (Zg-Nz/2).^2 15^2; rho(inner) 1.5; % 运行 SC 演化 2000 步 % 统计流量 rl mean(rho(inner rho0.4)); rv mean(rho(~inner rho0.4)); fprintf(密度比: %.2f\n, rl/rv);标定时要确认平衡后的密度不是初始给的1.5和0.8而是SC状态方程自己收敛出的值。密度比一旦低于目标档位加大G即可不用动其他参数。但G绝对值超过0.2后力尖峰会变得不可收敛此时应该换用更陡的伪势形如ψ ρ * exp(-ρ/ρ0)的变体工程上也比较常见。5.2 松弛时间τ与粘度、渗透率的换算τ决定流体运动粘度ν cs²(τ - 0.5)其中cs²1/3。多孔介质模拟对τ非常敏感τ太大边界层变厚微小孔隙内的流动被过度阻尼τ太接近0.5数值振荡增大。工程区间在0.7到1.5之间我自己一般固定在τ1.0起步验证通过后再改。用已知渗透率的解析阵列比如一排直圆柱组成的规则阵列标定模型渗透率时公式为K ν · Q · L / (Δt · ΔP)其中Q是体积流量、L是渗流长度。如果标定结果比解析值偏差超过5%先检查入口/出口缓冲区和力项速度修正这两个因素占误差来源大头。5.3 壁面吸附力Gads与接触角的关系接触角由固壁对两相的相对亲和性决定。Shan-Chen里常见的做法是在壁面格点施加额外吸附力Gads壁面处流体所受额外力与Gads符号相关Gads为正时壁面吸引流体中密度较高的部分接触角变小表现为亲液。标定接触角的标准手段是在一个平板通道内放一个静态液滴平衡后量三相接触点上的切向夹角。% 平板接触角标定z方向上下壁面 Gads 0.06; % 正亲湿相负憎湿相 % 在作用力函数里固体mask为1处额外加 for q 1:19 % 如果邻居是固体则把壁面密度当作常数rho_wall rho_wall 1.0; psi_w 1 - exp(-rho_wall); Fx Fx Gads * psi .* psi_w .* w(q) * c(q,1); end投影接触角与Gads的关系不是线性的通常要做三到四个不同Gads的标定实验拟合Gads-cosθ关系曲线。这样你在指定目标接触角后可以直接反查该用哪个Gads而不是盲猜之后反复跑全尺寸模拟。固体壁面的“刚度”也要注意如果壁面伪势太小界面临近壁面时会提前失稳出现不该有的附着层。5.4 时间步长与收敛判定SC模型的时间步长在LBM里是隐含的每个时间步对应物理时间Δt多用无量纲毛管数Ca来控制。多孔介质驱替里Ca μu/σ设目标Ca1e-4到1e-3反算驱替速度或者流量。运行时判断收敛不是看某个参数而是观察出口流量随时间是否进入平台期。这里给出常见的参数参考表表格里的量纲均为格子单位参数符号常用范围典型初值相互作用强度G-0.20 ~ -0.10-0.15平均密度ρ00.8 ~ 2.01.0松弛时间τ0.7 ~ 1.51.0壁面吸附力Gads-0.20 ~ 0.200.06伪势截断系数无1.0 ~ 2.01.0这些参数不是独立调节的。动了G之后密度比变化界面厚度也跟着变接触角标定会失准所以业界经验是先固定G调到目标密度比再固定τ调到目标粘度最后调Gads配接触角。中途换G前面标定全部重来。6. 验证与排错用Laplace标定法确认三维Shan-Chen实现可靠6.1 两球法测表面张力的最小代码实现SC之后第一件事不是跑多孔而是先在空旷网格里验证界面张力是否满足Laplace定律。做法是在网格中心放一个球形液滴等系统平衡后测液滴内外压差然后由ΔP 2σ/R反算σ。如果σ随R变化明显偏离常数说明力项离散有问题。% 静态液滴验证 R_list [12, 16, 20]; sigma_list zeros(3,1); for i 1:3 rho ones(Nx,Ny,Nz,single); R R_list(i); [Xg,Yg,Zg] ndgrid(1:Nx,1:Ny,1:Nz); inner (Xg-Nx/2).^2 (Yg-Ny/2).^2 (Zg-Nz/2).^2 R^2; rho(inner) 1.5; % 液滴初始密度 run_evolution(2000); % 取两相区域的体积平均压力 p_l mean(rho(inner) * rho0 / 3); p_v mean(rho(~inner) * rho0 / 3); sigma_list(i) (p_l - p_v) * R / 2; end if std(sigma_list) / mean(sigma_list) 0.05 fprintf(Laplace 验证通过σ≈%.4f\n, mean(sigma_list)); else warning(σ波动过大检查速度集和力项权重); end这个测试跑不通后面多孔模拟的毛管力全不可信。常见失败点是压力统计取了包含界面过渡层的格点统计区域要离界面至少两个格子否则压差偏大。6.2 三个高频排错点第一是质量守恒异常。每500步输出一次总质量sum(rho)若随时间单调增长或衰减超过0.5%要么是边界泄漏要么是Guo外力项在固体附近处理错误。第二是界面“碎成雾状”的抖动典型原因是G过大或τ小于0.7拉回参数区间即可如果是G值没问题但界面仍然抖动就检查伪势ψ是否在某些格点出现负值这种情况多半是密度初始扰动过猛。第三是驱替完全没有进展先看入口缓冲区是否被固体堵住再看连通域是否有中心孤岛导致压差传不进去。6.3 三维界面观察的一个实用技巧最后给一个我常用的MATLAB可视化检查法——不等时间步就调参数而是把驱替界面保存成isosurface数据若干步后再画。把f和rho的view只保留关键等值面会让内存占用小很多这一步用MATLAB自带的isosurface和patch画组相界面位置能直接看出前沿是否干净、是否有孤立液滴被拖后。% 用等值面提取两相界面 lv 0.5; % 密度等值面阈值两相密度中间值附近 [fv, fvc] isosurface(Xg, Yg, Zg, rho, lv); p patch(Faces, fv, Vertices, fvc); p.FaceColor interp; p.EdgeColor none; view(3); camlight; xlabel(x); ylabel(y); zlabel(z);这里核心技巧是不要直接画三维体数据体素太大MATLAB画图会卡到没法交互只用isosurface提取界面轮廓。界面形态从光滑曲面变成锯齿状碎片说明当前参数接近数值失稳边界及时降G或升高ρ0能救活很多本要重跑的长任务。参数没问题时这个界面应该会把孔隙通道逐步切开而不是在孔喉处留一条细长的灰色拖尾。本文还有配套的精品资源点击获取