简介这份资源聚焦岩土工程中基于LBM-DEM耦合方法的多相流固耦合数值模拟面向具备流体力学与固体动力学基础的研究生、科研人员及工程技术人员用于解决泥石流、颗粒沉降等复杂流固耦合问题。压缩包内包含1个PDF文档大小约843KB系统展示了论文复现思路、完整MATLAB代码及详细注释涵盖参数设置、LBM初始化、主循环、反弹边界与辅助函数等模块。内容不仅包括三个经典流固耦合案例的验证过程还深入分析了时间步长协调、多相界面交互等关键技术挑战揭示宏观与细观尺度上的流固相互作用机制。目前已有139人学习下载适合需要掌握LBM-DEM耦合原理、快速上手数值实现并开展二次开发的读者。该PDF将理论推导、代码实现与工程案例相结合可为相关方向的研究提供直接参考。1. LBM-DEM为什么在岩土流固耦合里比CFD-DEM更合适泥石流启动那几秒钟颗粒从静止渗流到突然失稳流体从孔隙里挤出来又裹着颗粒往下冲这个过程中最关键的细节都发生在颗粒尺度颗粒周围的绕流、颗粒间的碰撞、孔隙压力的瞬态变化。传统CFD-DEM方法把流体场离散成粗网格每个网格里可能包含几十个颗粒只能得到平均化的孔隙率和阻力颗粒表面附近的流动细节完全丢失。而格子玻尔兹曼方法LBM本身就是介观粒子模型网格可以细到和颗粒直径同一个量级颗粒边界在格子点上直接以浸入边界形式存在不需要反复重构体网格。离散元DEM这边由Cundall和Strack提出把岩土材料看成颗粒集合体颗粒间的接触力、摩擦、滚动都可以显式表达。两者耦合在一起正好补上了CFD-DEM在细观尺度上的短板。本文的落点是论文《Study on the multiphase fluid-solid interaction in granular materials based on an LBM-DEM coupled method》的复现骨架。论文用开源代码PalabosLBM和YadeDEM做了多相流固耦合模块并验证了瞬态流和泥石流两类岩土场景。我这里把它的核心逻辑用MATLAB重写成一个可以运行的简化版本覆盖D2Q9模型、反弹边界、拖曳力耦合和多相自由表面。适合有流体力学和固体力学基础的研究生、工程师目的是让你在动手写Palabos-Yade之前先用MATLAB把耦合流程跑通。2. D2Q9与LBM-DEM耦合机制从BGK碰撞到浸入边界2.1 D2Q9模型的九个方向与BGK近似LBM的核心不是直接求解纳维-斯托克斯方程而是追踪粒子分布函数 (f_i(x,t))表示在位置x处沿第i个离散方向运动的粒子数。D2Q9模型在二维平面内定义了9个离散速度方向中心静止方向0上下左右四个轴向方向以及四个对角方向。对应的权重依次为 (4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36)。权重的作用是保证平衡态分布函数在低马赫数下恢复出正确的宏观流体方程。分布函数的演化分为两步碰撞和流动。碰撞步用BGK单松弛近似把分布函数推向平衡态松弛时间 (\tau) 直接决定了流体的运动粘度关系式为[ \nu \frac{2\tau - 1}{6} ]这个公式从格子Boltzmann方程的多尺度展开得到(\tau 0.5) 对应零粘度实际计算中取值通常落在0.5到1.0之间。(\tau) 越接近0.5数值稳定性越差所以代码里用0.6或0.8是安全的选择。平衡态分布函数采用Maxwell-Boltzmann分布的离散形式[ f_i^{eq} \omega_i \rho \left(1 3(\mathbf{c}_i \cdot \mathbf{u}) \frac{9}{2}(\mathbf{c}_i \cdot \mathbf{u})^2 - \frac{3}{2}\mathbf{u} \cdot \mathbf{u}\right) ]宏观密度和速度通过矩计算得到(\rho \sum_i f_i)(\mathbf{u} \frac{1}{\rho}\sum_i f_i \mathbf{c}_i)。这段逻辑在MATLAB里对应compute_macroscopic和compute_equilibrium两个函数。2.2 流固信息传递的三种耦合方式LBM和DEM两个求解器如何交换信息是整个方法的要害。流场需要感受颗粒的存在颗粒需要受到流体的力和力矩。耦合方式基本思想精度实现难度适用场景拖曳力点耦合把颗粒简化为质点用经验阻力公式计算流体力低低颗粒尺寸远小于网格稀相流修正动量交换法在颗粒边界格子点上统计动量变化中中颗粒占据多个格子中等浓度浸入边界法IBM在碰撞算子中引入体积分数修正颗粒表面作为虚拟边界高高颗粒与流体强耦合密相流Noble和Torczynski最早把浸入边界法引入LBM-DEM思路是把固体体积分数 (\epsilon_s) 加进碰撞项颗粒内部的流体分布函数被强制向颗粒速度方向松弛。这样做的好处是颗粒不再是孤立的阻力点而是有空间占据的真实体。论文中Palabos和Yade的耦合模块采用的就是这类方法。MATLAB简化版里用的是第一类拖曳力耦合外加一步速度强制在颗粒占据的格子点上直接令流体速度等于颗粒速度然后重新计算平衡态分布函数。这样做在颗粒半径和网格尺度相当时会有一定的体积效应误差但作为框架演示足够清晰。2.3 多相流的自由表面简化处理真正的多相流固耦合需要同时追踪气相和液相界面常见做法是VOF或level set计算代价高。论文里针对泥石流场景简化处理为气液两相加上固体颗粒的三相问题。MATLAB代码里用一个分层密度函数模拟初始水位function rho handle_multiphase(rho, water_level, width, rho_air, rho_water) [nx, ny] size(rho); for y 1:ny if y water_level - width/2 rho(:,y) rho_air; elseif y water_level width/2 rho(:,y) rho_water; else alpha (y - (water_level - width/2)) / width; rho(:,y) rho_air alpha*(rho_water - rho_air); end end end这个函数按y坐标把计算域分成空气区、过渡区和水分区。过渡区宽度width控制界面的平滑程度width过小会导致密度突变在单松弛LBM中容易在界面处产生数值震荡过大则界面被抹得太宽失去多相意义。常见做法是把过渡区设为3到5个格子宽度。密度比也要注意。水和空气的真实密度比接近1000:1标准BGK-LBM在这个比值下会发散。MATLAB简化版里取0.1和1.0密度比只有10这是单松弛模型的稳定边界。要模拟真实水气比论文中的Palabos端采用的是多松弛时间MRT或者级数展开的界面力模型这也是LBM落地多相流时绕不开的进阶工作。3. MATLAB实现LBM-DEM主循环碰撞、流动、边界与颗粒更新3.1 参数初始化与格子单位LBM计算都在无量纲的格子单位下进行物理量要通过特征尺度换算。比如运动粘度的格子单位值 (\nu_{lattice}) 与物理值的关系是 (\nu_{physical} \nu_{lattice} \cdot \Delta x^2 / \Delta t)。这意味着网格间距和时间步长的选取直接决定了模拟对应多大尺度的问题。%% 参数设置 % LBM参数 nx 100; ny 100; % 网格尺寸 tau 0.6; % 松弛时间决定运动粘度 rho0 1.0; % 初始流体密度 u0 0.0; v0 0.0; % 初始流体速度 nu (2*tau - 1)/6; % 格子单位运动粘度 % DEM参数 particle_radius 5; % 颗粒半径格子单位 particle_density 2.0;% 颗粒密度 g 0.001; % 重力加速度格子单位 dt 1.0; % 时间步长这里particle_radius 5表示颗粒直径占10个格子这是浸入边界类方法能分辨颗粒形状的下限。如果颗粒半径只有1到2个格子颗粒在格子场里基本就是一个点必须改用拖曳力耦合。g 0.001看似很小但在格子单位下这相当于每时间步速度增加0.001个格子/步500步后颗粒速度也才0.5保证了低马赫数条件。如果直接用物理量纲的重力9.8LBM立刻发散。3.2 主循环的五个阶段主循环每个时间步顺序执行碰撞、流动、边界处理、颗粒更新、耦合反馈。其中流场的碰撞和流动是LBM的本体颗粒更新和耦合是DEM与LBM的连接点。for iter 1:max_iter % 1. LBM步骤 - 碰撞 [rho, ux, uy] compute_macroscopic(f); feq compute_equilibrium(rho, ux, uy, w, cx, cy); f f - (f - feq)/tau; % 2. LBM步骤 - 流动 f stream(f, cx, cy); % 3. 边界条件 - 反弹边界 f apply_bounce_back(f); % 4. DEM步骤 - 更新颗粒 [particle_x, particle_y, particle_vx, particle_vy] ... update_particle(particle_x, particle_y, particle_vx, particle_vy, ... rho, ux, uy, particle_radius, particle_density, g, dt); % 5. 耦合步骤 - 颗粒速度传回流体 f apply_coupling(f, particle_x, particle_y, particle_radius, ... particle_vx, particle_vy, rho, ux, uy); % 可视化 if mod(iter, 10) 0 visualize(rho, particle_x, particle_y, particle_radius, iter); end end碰撞步的逻辑是f f - (f - feq)/tau含义是把当前分布函数向平衡态推进1/tau的比例。流动步stream把每个方向的分布函数沿离散速度方向移动到相邻格子。这里的mod索引实现了周期性边界粒子从右边界流出就从左边界流入。但岩土流固耦合通常需要封闭边界所以在后面的步骤中会用反弹边界覆盖掉周期边界的行为。值得注意的一个坑是compute_macroscopic函数内部直接引用了变量cx和cy而这两个变量在函数定义中没有被作为参数传入。如果函数保存在独立脚本文件里MATLAB会报Unrecognized function or variable cx。有两个解决办法一是把速度数组直接写死在函数内部二是把cx、cy放进结构体params所有函数都接收params作为参数。第二种方式更干净也方便批量调整参数。3.3 反弹边界的索引映射反弹边界是LBM里处理固壁最简单的方式粒子撞到壁面后沿原路返回。落实在代码上就是把对应方向的分布函数对调。标准反弹映射关系为方向1对32对45对76对8。function f apply_bounce_back(f) % 下边界y1处方向4、7、8 反弹为 2、5、6 f(2:end-1, 1, [4,7,8]) f(2:end-1, 1, [2,5,6]); % 上边界yny处方向2、5、6 反弹为 4、7、8 f(2:end-1, end, [2,5,6]) f(2:end-1, end, [4,7,8]); % 左边界x1处方向3、6、7 反弹为 1、5、8 f(1, 2:end-1, [3,6,7]) f(1, 2:end-1, [1,5,8]); % 右边界xnx处方向1、5、8 反弹为 3、6、7 f(end, 2:end-1, [1,5,8]) f(end, 2:end-1, [3,6,7]); end注意这里的索引写法f(2:end-1, 1, [4,7,8])只处理了上下边界的中段四个角落没有做处理。角落同时属于两条边需要额外定义角点规则。对于颗粒沉降这类中心区域流动的问题角落误差对结果影响有限但如果模拟管道流或者有障碍物的流场角落处理不能省略。3.4 颗粒受力的拖曳力模型update_particle函数实现了颗粒在流场中的受力更新。流体对颗粒的作用力用修正后的拖曳力公式计算阻力系数依赖雷诺数% 计算拖曳力系数 Cd 0.4; Re 2*radius * sqrt((vx-ux_f)^2 (vy-uy_f)^2) / 0.1; if Re 0 Cd 24/Re 6/(1sqrt(Re)) 0.4; end % 拖曳力 Fd_x 0.5*Cd*rho_f*pi*radius^2 * (ux_f-vx) * abs(ux_f-vx); Fd_y 0.5*Cd*rho_f*pi*radius^2 * (uy_f-vy) * abs(uy_f-vy); % 浮力修正后的重力 Fg (density - rho_f) * pi * radius^2 * g; % 更新速度和位置 vx_new vx Fd_x * dt / (density*pi*radius^2); vy_new vy (Fd_y Fg) * dt / (density*pi*radius^2);雷诺数计算里除以0.1用的是格子单位运动粘度的假定值。这里代码直接硬编码了0.1如果前面tau改了这个值必须同步改否则阻力系数和实际流场粘度不一致。一个更规范的做法是把nu作为参数传入这个函数。阻力公式中的abs(ux_f-vx)是为了保证阻力的方向性颗粒相对流体运动方向不同阻力方向也不同。Fg用的是颗粒密度减流体密度再乘体积和重力相当于自动考虑了浮力不需要再额外写浮力项。这是颗粒沉降模拟的标准写法。4. 从单颗粒到多颗粒接触力、参数标定与时间步协调4.1 多颗粒系统的结构体组织与接触力单颗粒版本只能验证耦合算法流程真正要模拟泥石流或瞬态流需要几十上百个颗粒。输入代码里给出了一个5颗粒版本数据结构从独立变量改成结构体数组每个颗粒带位置、速度、半径、密度。规模扩大后颗粒间的接触力就变得和流固耦合同等重要。DEM部分的核心是接触力学模型。颗粒之间发生重叠时法向接触力用线性弹簧模型% 颗粒i和颗粒j之间的距离 dx particles(j).x - particles(i).x; dy particles(j).y - particles(i).y; dist sqrt(dx^2 dy^2); % 重叠量 overlap particles(i).radius particles(j).radius - dist; if overlap 0 % 法向接触力 Fn kn * overlap; % 切向接触力简化的库仑摩擦 Fs kt * tangential_displacement; if Fs mu * Fn Fs mu * Fn; end endkn是法向接触刚度kt是切向刚度mu是摩擦系数。刚度的选取直接决定接触时间尺度。根据弹性接触理论两个颗粒碰撞的接触时间约为 (\pi \sqrt{m/k_n})。这个量必须和LBM的时间步长匹配否则颗粒碰撞产生的力波无法被流场正确响应。参数物理含义取值范围建议对模拟的影响kn法向接触刚度使最大重叠量小于颗粒半径的1%过大导致接触时间过短过小颗粒发生明显穿透kt切向接触刚度通常取0.2~0.5倍kn影响颗粒旋转和摩擦行为mu库仑摩擦系数砂土取0.3~0.6黏土取0.1~0.3影响堆积角、休止角g重力加速度格子单位下需保证最大流速小于0.1决定颗粒沉降速度和流场稳定性多颗粒版本在LBM循环中用update_particles代替单颗粒的update_particle内部嵌套两层循环检测每一对颗粒的接触。时间复杂度是O(N^2)颗粒数超过几百个后计算量急剧上升。Yade的工程实现里用的是网格邻居搜索和Verlet列表MATLAB版本做好O(N^2)的验证就够了真正的百万颗粒规模还是要靠C。4.2 时间步长协调策略LBM-DEM耦合最常见的不稳定来源是时间步长不匹配。LBM的时间步被格子粘度和CFL条件限制要求格子速度 (u_{max} 0.1) 左右DEM的时间步受接触刚度和颗粒质量限制要求 (\Delta t_{DEM} \pi\sqrt{m/k_n})。两者的约束条件完全不同直接共用一个时间步长往往导致要么流场发散要么颗粒穿透。常见的工程做法是子循环LBM走一步的时间段内DEM走若干子步。比如LBM步长 (dt_{LBM} 1)DEM子步 (dt_{DEM} 0.1)耦合界面每10个DEM子步同步一次流场数据。这样做的好处是流固耦合的信息交换频率保持不变但颗粒接触的稳定性得到保证。需要注意避免的是在两个求解器之间来回插值过度每个子步都重新读取流体速度会导致数值耗散累积。另一个时间尺度问题是颗粒曳力响应时间。颗粒在流场中达到力平衡的特征时间约为 (t_{relax} \frac{\rho_p d_p^2}{18\mu_f})。如果这个时间远大于一个LBM时间步说明颗粒对流场的响应是缓慢的可以降低耦合频率每5到10个LBM步同步一次。如果这个时间小于时间步长说明颗粒运动太快需要减小时间步或者加密网格。论文中泥石流模拟的尺度范围从稳定渗流阶段到启动阶段再到流固混合物流动阶段三个阶段的时间尺度相差很大。稳定渗流时流体速度极低颗粒几乎不动如果全程用同一个时间步长要么前期计算浪费要么后期发散。实际工程中会在不同阶段切换耦合频率甚至切换网格分辨率。4.3 多相界面的参数调试多颗粒版本的代码里多了一个handle_multiphase函数把流场初始化为上空气下水的分层结构。这个初始化的意义是模拟水库水位骤降或降雨入渗后形成的非饱和区域。interface_width和water_level两个参数控制初始界面的位置和过渡带厚度。调试时重点关注三个现象。第一界面处密度突变会产生初始压力震荡表现为密度场出现波纹状图案这时应增大interface_width而不是减小tau。第二颗粒穿过界面时阻力系数的雷诺数计算要确保使用颗粒所在位置的局部流体参数而不是全局平均值。第三如果空气区密度设置过低会导致空气区速度异常增大因为低压区的微小动量就会产生高速。取rho_air 0.1而rho_water 1.0时空气区速度大约是水区的3倍这在LBM中是可以接受的。下面给出多颗粒版本初始化时的一个关键片段% 多相流参数 rho_air 0.1; % 空气密度 rho_water 1.0; % 水密度 interface_width 3; % 界面过渡宽度格子数 water_level ny/2; % 初始水位 % 颗粒初始位置避开界面按列分布避免初始接触重叠 for i 1:num_particles particles(i).x 20 (nx-40)*rand(); particles(i).y ny - 20*i; end颗粒初始化时用ny - 20*i将颗粒从顶部向下等间距排列。这样做避免了随机分布导致初始重叠省去了DEM接触检测在第一个时间步就产生巨大排斥力的麻烦。rand()只用于水平位置扰动保证颗粒不会全部落在同一垂直线上。5. 验证、后处理与Palabos-Yade落地路线5.1 用颗粒终速公式验证耦合正确性LBM-DEM代码跑通之后第一个动作不是调参优化而是验证物理量是否合理。单颗粒沉降问题有解析解颗粒在流体中达到平衡时重力、浮力和拖曳力满足力平衡终速为[ v_t \sqrt{\frac{4 g d_p (\rho_p - \rho_f)}{3 C_d \rho_f}} ]其中 (d_p 2 \times particle_radius)(C_d) 通过雷诺数迭代计算。可以在MATLAB里做一个独立脚本把主循环输出的颗粒速度与这个解析值对比。误差在5%到10%是可接受的因为格子单位下离散化和边界处理会带来微小偏差。% 后处理脚本片段读取最后100步的颗粒速度并求平均 load(particle_trajectory.mat); terminal_velocity mean(vy(end-100:end)); Cd 24/Re 6/(1sqrt(Re)) 0.4; vt_theory sqrt(4*g*(2*particle_radius)*(particle_density - rho0)/(3*Cd*rho0)); error_pct abs(terminal_velocity - vt_theory) / vt_theory * 100; fprintf(仿真终速: %.4f, 理论终速: %.4f, 误差: %.2f%%\n, ... terminal_velocity, vt_theory, error_pct);如果误差偏大优先检查雷诺数计算中的运动粘度是否与tau一致。代码里update_particle函数中雷诺数计算硬编码了0.1这个粘度值这是最容易出错的地方。5.2 密度场可视化与数据导出主循环里的visualize函数用imagesc绘制密度场rectangle叠加圆形表示颗粒位置。这个方法在调试时直观高效但要输出论文级别的图还需要改进用surf加光照显示界面形状或者用contourf画等密度线。另外建议每个时间步把颗粒位置和流场密度保存到.mat文件后处理阶段再做动画和定量分析。使用save导出数据时要注意主循环运行中直接调用save会产生大量IO开销。常见做法是先预分配数组在循环内只做内存赋值循环结束后统一保存trajectory zeros(max_iter, 2); for iter 1:max_iter % ... 主循环 ... trajectory(iter, :) [particle_x, particle_y]; end save(particle_trajectory.mat, trajectory);5.3 从MATLAB原型到Palabos与YadeMATLAB代码的定位是原型验证它跑通的是耦合逻辑和参数量级。论文中的实际算例用的是Palabos和Yade两个开源框架。Palabos负责LBM部分采用C实现支持多松弛时间和多相流扩展Yade是DEM引擎支持复杂颗粒形状和接触本构模型。两者的数据交换通过socket或者共享内存实现每个耦合步交换颗粒位置、速度和流体作用力。从MATLAB向Palabos-Yade迁移时有三个关键差异需要提前准备。第一Palabos的网格是所有进程共享的分布式网格颗粒覆盖多个进程时需要MPI层面的通信。第二Yade的接触模型比MATLAB版丰富得多包括Hertz-Mindlin、JKR黏聚接触等参数含义和MATLAB版不同。第三真实的多相流模型需要引入表面张力项其离散形式在LBM中表现为额外的力项这个在MATLAB简化版里完全被省略了。如果只是做小规模验证第5.1节的解析解验证方法可以直接复用到Palabos输出上。建议先在MATLAB里把沉降问题调到误差5%以内再迁移到Palabos和Yade这样可以避免两个大型代码库同时排错的困境。对于需要进一步加速的场景可以将stream循环写成mex文件或者用parfor并行化这能获得数倍提速同时也让结果与实验结果对比时更具说服力。本文还有配套的精品资源点击获取