简介基于LBM-DEM耦合方法的多相流固耦合模拟资料聚焦岩土工程中的复杂流固相互作用问题适合具备流体力学与固体动力学基础的硕博研究生及从事泥石流、颗粒沉降等数值模拟的工程师。内容以论文复现为主线提供可运行的MATLAB代码涵盖参数设置、LBM初始化、宏观量计算、平衡分布函数、碰撞与流步、反弹边界及颗粒位置更新等核心模块并在代码注释中解释每一步的物理含义。全文围绕瞬态流与泥石流等经典案例展开同时探讨了时间步长协调、多相界面交互等关键难点有助于理解细观尺度的流固耦合机制。全部资料为一个843KB的PDF文件结构紧凑既有算法框架说明也有代码逐段注释便于读者快速上手和二次开发。该资源目前已有139人浏览学习对希望掌握LBM-DEM耦合原理并用于工程验证的读者具有直接的参考价值。1. 岩土工程里的LBM-DEM耦合为什么值得自己写一遍岩土工程里的管涌、流土、降雨入渗和泥浆渗透本质上都是多相流固耦合问题水、气两相在孔隙里流动土颗粒被冲刷、搬运、重新排列形成新的孔道。传统连续介质方法把土骨架固定死用达西定律算渗透一旦颗粒开始运动渗透系数、孔隙率和接触网络全部在变方程本身就不成立了。LBM-DEM耦合是把两套思路拼在一起LBM格子玻尔兹曼方法在孔隙尺度解流动复杂孔隙边界用反弹格式直接处理绕开有限元里反复网格重构的麻烦DEM离散元负责颗粒接触和运动。两者通过动量交换把流体力传给颗粒、把颗粒位移反馈给流场就能模拟水流搬运颗粒、颗粒堵住孔道、孔道改变流场的完整闭环。这套方法适合数值分析方向的岩土研究生、做颗粒材料模拟的工程师以及想把论文算法落地成可运行代码的人。MATLAB 在这里是顺手的原型平台LBM 的局部更新特性让矩阵化实现很紧凑不需要优化工具箱基础版就能跑通。下面按理论映射、最小代码、参数与验证、落地技巧四层拆开讲。2. LBM-DEM耦合的理论骨架格子玻尔兹曼、离散元与力传递的三层映射先把三层各自干的事理清LBM 解流体相存的是分布函数而不是压力速度场DEM 解离散颗粒相存的是颗粒位置、速度与接触力两者之间的公共接口是颗粒表面处流体的拖曳力和颗粒边界对流场的反射作用。所有量在实现阶段都先用格子单位无量纲跑通之后再统一换算到物理单位这一步是后面所有参数讨论的前提。2.1 D2Q9格子与BGK碰撞流场状态怎么存、怎么更新LBM 不解 N-S 方程本身而是在每个格点存 9 个方向的分布函数 f_k(x,t)。每一步只有两个操作碰撞把 f_k 向平衡态分布 f_k^eq 松弛松弛快慢由无量纲松弛时间 tau 控制迁移把每个 f_k 沿其速度方向 e_k 平移一个格点。宏观密度和动量是 f_k 的零阶矩和一阶矩ρ Σf_kρu Σf_k e_k。二维最常用 D2Q99 个方向的权重和速度矢量是固定的k方向 (ex, ey)权重 w_k1(0, 0)4/92(1, 0)1/93(0, 1)1/94(-1, 0)1/95(0, -1)1/96(1, 1)1/367(-1, 1)1/368(-1, -1)1/369(1, -1)1/36在 MATLAB 里平衡态分布就是一个向量化函数这是整个耦合代码里最基础的积木function feq equilibrium(rho, ux, uy, w, ex, ey) % D2Q9 平衡态分布cs2 1/3 时系数的简化形式 feq zeros([size(rho), 9]); for k 1:9 eu ex(k)*ux ey(k)*uy; % e_k · u usq ux.^2 uy.^2; % u · u feq(:,:,k) w(k) .* rho .* ... (1 3*eu 4.5*eu.^2 - 1.5*usq); end end这里的核心参数是 tau。运动粘度在格子单位下等于 ν (tau - 0.5)/3所以 tau 必须大于 0.5 才是正粘度。tau 越接近 0.5数值粘度越低流场越容易振荡tau 超过 1.2 后数值耗散又过大孔隙里的低速流动会被抹平。岩土渗流这种低雷诺数场景tau 取 0.71.0 最稳。2.2 DEM颗粒接触模型法向弹簧阻尼与切向库仑摩擦DEM 把每个土颗粒当作刚性圆盘二维或球三维颗粒之间用弹簧-阻尼器接触模型。法向力只压不拉切向力用弹簧加库仑摩擦截断。这是 LIGGGHTS、EDEM 这类软件里线性接触模型的公共内核MATLAB 里写成一个独立函数最合适function [Fn, Ft] demContact(delta_n, vn, delta_t, vt, kn, cn, kt, ct, mu) % 法向max(0, ...) 保证只压不拉 % vn 取“分离速度”因此阻尼项带负号 Fn max(0, kn * delta_n - cn * vn); % 切向弹性 阻尼再用库仑摩擦截断 Fte kt * delta_t - ct * vt; Ft sign(Fte) .* min(abs(Fte), mu * abs(Fn)); end参数之间有硬约束法向刚度 k_n 直接决定允许的最大重叠量 δ_n F/k_n工程经验是重叠量控制在颗粒直径的 1% 以内k_t 一般取 (2/7 ~ 1/2) k_n阻尼系数 c_n 通常由恢复系数 e 反推c_n -2 ln(e)√(m_eff k_n / (ln²e π²))。稳定性约束是 DEM 时间步必须小于接触振荡周期的 2√(m/k_n)实际取该值的 0.2 倍以内。2.3 多相流的Shan-Chen伪势模型与流固耦合的力接口多相部分最常见的做法是 Shan-Chen 伪势模型在碰撞前额外加一个体积力力的大小由局部伪势梯度决定。伪势取 ψ(ρ) ρ0(1 - e^(-ρ/ρ0))作用强度 G 是负值驱动同相颗粒聚集界面张力和两相密度比都由 G 与状态方程共同控制。G 的绝对值稍大于 1 就开始分相实际常用区间见第 4 章。流固耦合的力接口用半程反弹加动量交换检测流体节点与固体节点之间的链接把指向固体的分布函数反弹回反方向固体颗粒一次收到 2·e_k·f_k 的动量。所有边界链接的动量累加就是颗粒受到的流体力。颗粒移动一步后重新标记固体节点孔隙几何随颗粒位移实时更新这就是复杂流固相互作用里几何不断变化的那一层。另一种常见接法是浸没边界法IBM力通过插值扩散到流场网格不需要贴合颗粒表面但对岩土里大量不可渗透颗粒反弹格式更直观孔隙结构直接由一张标记矩阵表达。3. MATLAB实现LBM-DEM多相流固耦合最小可运行代码与逐段说明下面这套代码能跑通流体搬运颗粒、颗粒阻塞孔道的完整闭环。它基于二维 D2Q9不依赖 Simulink 和优化工具箱MATLAB R2020b 之后的版本直接可以运行。代码刻意写成教学骨架每个数组都保持完整形状方便逐步打印检查。3.1 数据组织颗粒、流场与固体标志位的矩阵化安排初始化把所有状态量安排成三个部分分布函数 f 是 Nx×Ny×9 的三维数组第三维是方向这是后续所有矩阵化操作的支点颗粒状态用一组行向量存位置、速度和力固体标志位 solidMat 用整数矩阵0 表示流体正整数表示颗粒编号反弹时能直接定位到具体颗粒。% 网格与流体参数格子单位 Nx 300; Ny 150; % 流场网格数 tau 0.8; omega 1/tau; cs2 1/3; % tau 控制粘度 % D2Q9 方向与权重 ex [0, 1, 0, -1, 0, 1, -1, -1, 1]; ey [0, 0, 1, 0, -1, 1, 1, -1, -1]; w [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; opp [1, 4, 5, 2, 3, 8, 9, 6, 7]; % 反弹方向索引 % 分布函数初始化为平衡态rho1, u0 f zeros(Nx, Ny, 9); for k 1:9 f(:,:,k) w(k); end rho ones(Nx, Ny); ux zeros(Nx, Ny); uy zeros(Nx, Ny); % 颗粒数据 Np 12; R 5; % 颗粒数、半径格点 xp 10 (Nx-20)*rand(1,Np); % 随机布点避开边界 yp 10 (Ny-20)*rand(1,Np); vxp zeros(1,Np); vyp zeros(1,Np); Fpx zeros(1,Np); Fpy zeros(1,Np); rho_p 2.0; % 颗粒相对密度 % 固体标志位0 流体p 颗粒编号 solidMat zeros(Nx, Ny); [xx, yy] meshgrid(0.5:1:Nx-0.5, 0.5:1:Ny-0.5); for p 1:Np solidMat((xx-xp(p)).^2 (yy-yp(p)).^2 R^2) p; end初始化里最值得注意的一点固体节点上的分布函数初始也是 w_k没有特殊处理。流体被颗粒占住的位置在第一步迁移后会被清零见 3.2 第 6 步这套先统一初始化、运行中持续清零的写法比只初始化流体节点简单得多也方便调试时直接检查 solidMat 的形状。3.2 多相流LBM核心伪势力、碰撞、反弹与迁移的矩阵写法单步 LBM 的顺序是宏观量 - 伪势力 - 速度修正 - 碰撞 - 反弹 - 迁移。反弹必须放在迁移之前否则分布函数已经进入固体节点反弹就失去了物理意义动量交换力在反弹的同时累加省掉再扫一遍网格。% LBM 单步主循环内调用 % 1) 宏观量对第三维求和/加权 rho sum(f, 3); ux sum(f .* reshape(ex,1,1,9), 3) ./ rho; uy sum(f .* reshape(ey,1,1,9), 3) ./ rho; % 2) Shan-Chen 伪势力 rho0 1.0; Gsc -1.8; % Gsc 为负驱动分相 psi rho0 * (1 - exp(-rho / rho0)); Fsc_x zeros(Nx, Ny); Fsc_y zeros(Nx, Ny); for k 1:9 ps circshift(psi, [ex(k), ey(k)]); Fsc_x Fsc_x - Gsc * w(k) .* psi .* ps .* ex(k); Fsc_y Fsc_y - Gsc * w(k) .* psi .* ps .* ey(k); end % 3) 伪势力折算进平衡态速度Shan-Chen 原方案的做法 ux ux tau .* Fsc_x ./ rho; uy uy tau .* Fsc_y ./ rho; % 4) 碰撞向平衡态松弛 for k 1:9 eu ex(k)*ux ey(k)*uy; feq w(k) .* rho .* ... (1 3*eu 4.5*eu.^2 - 1.5*(ux.^2 uy.^2)); f(:,:,k) (1 - omega) * f(:,:,k) omega * feq; end % 5) 半程反弹 动量交换放在迁移之前 for k 1:9 sNbr circshift(solidMat, [-ex(k), -ey(k)]); % k 方向邻居 swap (solidMat 0) (sNbr 0); % 流体节点紧邻颗粒 tmp f(:,:,k); % 交换正反方向分布实现反射 f(:,:,k) tmp .* ~swap f(:,:,opp(k)) .* swap; f(:,:,opp(k)) f(:,:,opp(k)) .* ~swap tmp .* swap; % 动量交换固体收到 2*e_k*f_k静止壁面近似 [ii, jj] find(swap); for n 1:numel(ii) p sNbr(ii(n), jj(n)); % 归属颗粒编号 Fpx(p) Fpx(p) 2 * ex(k) * tmp(ii(n), jj(n)); Fpy(p) Fpy(p) 2 * ey(k) * tmp(ii(n), jj(n)); end end % 6) 迁移周期性边界下直接移位固体节点清零 for k 1:9 f(:,:,k) circshift(f(:,:,k), [ex(k), ey(k)]); f(:,:,k) f(:,:,k) .* (solidMat 0); end这段代码里三个细节直接决定结果对不对。第一伪势力用circshift(psi, [ex(k), ey(k)])取邻点伪势迁移也是同一套位移两处方向约定必须完全一致否则力和流动方向错位。第二第 5 步的find循环只遍历边界链接数量是颗粒周长的量级不会拖慢整体真正的性能热点在碰撞和迁移。第三第 6 步的迁移用了circshift它自带周期性边界如果要做上下壁面、左右压差的渠道流需要把上下边界替换成反弹壁面不能直接用这行。注意上面是静止壁面近似。颗粒在运动时反弹分布要减去壁面速度贡献tmp - 2*w(k)*rho.*(ex(k)*vxp(p)ey(k)*vyp(p))/cs2否则颗粒运动会往流场注入虚假动量算出来的沉降速度会偏大。3.3 DEM步进与固体掩膜更新从流体力到位移的完整闭环LBM 输出了作用在每个颗粒上的流体力 Fpx/FpyDEM 部分要做三件事叠加重力和浮力、算颗粒间接触力、用速度 Verlet 更新位移。主循环里每步开始要把 Fpx、Fpy 清零因为它们在上一步里已经消耗掉了。% DEM 一步净重力 接触力 速度 Verlet g_lb 5e-5; % 格子单位重力向下为正 mass rho_p * pi * R^2; % 二维单位厚度圆盘 % 净重力 (颗粒密度-流体密度)*体积*g静水压已经含在动量交换里 Fpy(:) Fpy(:) - (rho_p - 1) * pi * R^2 * g_lb; % 第一步半步速度 全量位移LBM 步长 dt1 vxp vxp 0.5 * Fpx / mass; vyp vyp 0.5 * Fpy / mass; xp xp vxp; yp yp vyp; % 颗粒两两接触法向线性接触切向见 2.2 节的函数 for p 1:Np for q p1:Np dx xp(q)-xp(p); dy yp(q)-yp(p); d sqrt(dx^2dy^2); n [dx, dy]/d; delta_n 2*R - d; % 法向重叠量 if delta_n 0 vn (vxp(q)-vxp(p))*n(1) (vyp(q)-vyp(p))*n(2); Fn max(0, kn*delta_n - cn*vn); Fpx(p) Fpx(p) - Fn*n(1); Fpy(p) Fpy(p) - Fn*n(2); Fpx(q) Fpx(q) Fn*n(1); Fpy(q) Fpy(q) Fn*n(2); end end end % 第二步完成速度更新加一点数值阻尼抑制颗粒抖振 vxp vxp 0.5 * Fpx / mass; vyp vyp 0.5 * Fpy / mass; vxp vxp * 0.999; vyp vyp * 0.999; % 用新位置重建固体掩膜 solidMat zeros(Nx, Ny); for p 1:Np solidMat((xx-xp(p)).^2 (yy-yp(p)).^2 R^2) p; end这段的物理关键是净重力表达式。动量交换已经在流体力里包含了压力项和粘性力项所以重力必须用浮容重即颗粒密度减流体密度再乘以体积否则颗粒即使静止也会因虚力缓慢下沉。接触部分的双重循环在 Np50 以下完全够用颗粒数上千时再改成 cell list 或邻接表。颗粒移开后暴露出来的原固体节点下一轮 LBM 会从邻居迁移进分布函数但为了不出现局部质量缺失正式计算时要在重建 solidMat 后把这些节点用周围流体的平均密度和速度重新初始化。主循环就是把 3.2 和 3.3 按顺序串起来每步先清零 Fpx/Fpy做 LBM 单步含反弹与动量交换再做 DEM 接触与运动更新最后重建 solidMat。跑 10000 步左右颗粒沉降和两相分布就能进入统计稳态。4. 复杂流固相互作用的参数设定与算例验证松弛时间、润湿性与颗粒刚度耦合模型能跑起来是一回事参数对不上是另一回事。多相 LBM 与 DEM 各有自己的数值限制合在一起时约束叠加。经验顺序是先固定流体参数再标定两相界面参数最后校核接触刚度与时间步。4.1 松弛时间、密度比与Shan-Chen作用强度的设定顺序先定 tau。岩土渗流多为低雷诺数流动tau 取 0.7 附近时数值粘度小流场对颗粒运动的响应快但需要配合较细的网格tau 取 1.0 时非常稳代价是渗透率测量结果偏保守。再定 Gsc标准指数型伪势下|Gsc| 稍大于 1 就开始分相实际常用 -1.2 到 -2.0。Gsc 绝对值越大界面张力越高、两相密度比越大但界面处伪速度spurious velocity也越剧烈颗粒表面附近的密度翼会污染流固耦合力。关于密度比要泼一盆冷水标准 Shan-Chen 指数伪势撑死做到几十比一的密度比直接做水-气约 1000:1会失稳。两相密度接近的体系比如水-非水相液体在孔隙里的驱替用标准模型最合适非要做水-气需要换 Carnahan-Starling 状态方程或用颜色梯度模型。润湿性通过额外加一项流-固伪势力实现系数符号决定接触角亲水取负、疏水取正这是复杂流固相互作用里最容易忘的一层。参数符号建议初值范围主要影响取值过当的后果松弛时间tau0.6 ~ 1.0粘度与流场稳定性小于 0.55 高频振荡大于 1.2 耗散过大伪势作用强度Gsc-1.2 ~ -2.0相分离强度、界面张力过小不分相过大界面伪速度失控伪势参考密度rho01.0两相密度比上限水-气场景必须换状态方程法向接触刚度k_n1e4 ~ 1e6颗粒重叠量、DEM 步长太小重叠不可忽略太大幅度受限切向刚度k_t(2/7 ~ 1/2)k_n摩擦角与堆积形态过大引发接触抖动颗粒半径R5 ~ 10 格点孔隙分辨率、渗透率标定小于 3 格点标定结果失真4.2 颗粒接触刚度与时间步长的匹配约束DEM 与 LBM 耦合时时间步的约束是叠加的LBM 要求 tau 0.5 且马赫数足够低DEM 要求步长远小于接触振荡周期。实际操作里先把 k_n 定下来估算单颗粒受到的净重力加流体拖曳力 F_est要求重叠量 δ_n F_est/k_n 小于颗粒直径的 1%由此反推 k_n 的下限。然后计算 DEM 临界步长 Δt_dem 2√(m_eff/k_n)取 0.2 倍再与 LBM 步长比较。如果 Δt_dem 远小于 LBM 步长两个选择一是把 k_n 调小到两者同量级二是在一个 LBM 步内做若干次 DEM 子循环。经验是纯教学代码用方案一更省事k_n 取 1e4 量级就能让重叠量在 5% 半径以内要模拟硬颗粒重叠量小于 1%必须用方案二。另一个常被忽略的点接触阻尼不能给太大c_n 过大会让颗粒像泡在糖浆里沉降速度系统性偏小用恢复系数 e 反推 c_n 是最稳的做法。4.3 两个经典验证算例单颗粒沉降与达西渗流标定验证是参数标定的前提。第一个算例是单颗粒沉降二维模型必须注意圆盘在纯蠕流里没有 Stokes 解要用 Oseen 修正后的圆筒阻力公式 F/L 4πμU / (0.5 - γ - ln(Re/8))γ 0.5772。用这个公式解出终端速度与数值结果对比。严格对照 Stokes 沉降必须上三维 D3Q19二维对比只能用阻力系数公式或者 Re 1 时的 Schiller-Naumann 修正。第二个算例更贴近岩土让单相流体穿过随机堆积的颗粒床测达西渗透率。稳态后在模型中部取一个窗口平均流速压力梯度用多相状态方程的压力差算而不是直接拿密度差换算% 稳态后测渗透率 p_field cs2*rho cs2*Gsc*psi.^2/2; % SC 状态方程压力 dpdx (mean(p_field(1,:)) - mean(p_field(Nx,:))) / Nx; u_darcy mean(mean(ux(Nx/3:2*Nx/3, Ny/4:3*Ny/4))); % 中部窗口平均 K_lb u_darcy * ((tau-0.5)/3) / dpdx; % 格子渗透率 % 换算到物理单位并对比 Kozeny-Carman 参考值 dx 1e-4; % 每个格点对应的物理尺寸米 K_phys K_lb * dx^2; K_kc d_p^2 * phi^3 / (180 * (1-phi)^2); fprintf(数值 K %.3eKozeny-Carman %.3e\n, K_phys, K_kc);注意 Kozeny-Carman 是从三维堆积推导的经验公式二维模型里只能当量级参考。真正的标定做法是把 dpdx 换成一个已知渗透率的标准多孔介质反推出 dx 和 tau 的组合然后用这个组合去算目标工况。窗口位置也要固定颗粒重排后孔隙率随高度变化窗口选在中部 1/3 高度处最稳定选在边界附近测出的渗透率会差一个量级。5. 岩土工程应用里的三个落地技巧向量化加速、渗透率标定与可视化输出5.1 用circshift和第三维广播代替三层循环初学版本最常见的性能杀手是按格点双重循环算碰撞。矩阵化的核心是让第三维9 个方向直接参与广播运算sum(f,3)替代for i, for j的两层聚合circshift替代迁移的逐格点搬运。碰撞还能进一步改写成矩阵运算把 f 重排成 Nx*Ny 行 9 列的矩阵平衡态计算一次性完成代码同样简洁。显存够用时把 f、rho、ux 用gpuArray包一层碰撞和迁移的代码一行都不用改规模在 500×500 网格以下能获得明显加速。5.2 从格子单位换算到物理单位渗透率标定的一条路径换算链只有三步先定 dx一个格点对应的物理长度再由物理粘度 ν_phys 和格子粘度 ν_lb (tau-0.5)/3 反推时间步 dt dx²·ν_lb/ν_phys最后任何格子渗透率 K_lb 都按 K_phys K_lb·dx² 转换。工程上更常用的是反向标定手头有一组室内渗透试验的 K_exp固定 dx调节 tau 代进数值模型跑一次达西算例直到 K_phys 与 K_exp 吻合这一步实际是把网格尺度和数值粘度同时吸收进了模型参数后续所有工况都用这一组换算系数。5.3 后处理孔隙水压力、两相分布与颗粒位移的同步可视化MATLAB 画图的后处理三段式两相分布用imagesc(rho)并叠加阈值等值线压力场用状态方程算完再contourf颗粒用viscircles叠加在图上。颗粒细小、数量多时只画抽样位置的 quiver全部矢量画出来图面会糊。每 N 步存一帧最后用VideoWriter合成动画比实时绘图省内存。测量任何宏观量时都记住一个原则取模型中部固定窗口的空间平均不要取全流域。颗粒重排、入口效应和出口回流都会污染统计值中部窗口是唯一能同时避开这三类干扰的位置。本文还有配套的精品资源点击获取