企业官网建设流程全解析

简介：基于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。

二维最常用 D2Q9，9 个方向的权重和速度矢量是固定的：

在 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.7～1.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]; % 反弹方向索引 % 分布函数初始化为平衡态：rho=1, u=0 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/Fpy，DEM 部分要做三件事：叠加重力和浮力、算颗粒间接触力、用速度 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 步长 dt=1） 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 = p+1:Np dx = xp(q)-xp(p); dy = yp(q)-yp(p); d = sqrt(dx^2+dy^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

这段的物理关键是净重力表达式。动量交换已经在流体力里包含了压力项和粘性力项，所以重力必须用浮容重，即颗粒密度减流体密度再乘以体积，否则颗粒即使静止也会因虚力缓慢下沉。接触部分的双重循环在 Np=50 以下完全够用，颗粒数上千时再改成 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 状态方程或用颜色梯度模型。润湿性通过额外加一项流-固伪势力实现，系数符号决定接触角：亲水取负、疏水取正，这是复杂流固相互作用里最容易忘的一层。

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 = %.3e，Kozeny-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合成动画，比实时绘图省内存。测量任何宏观量时都记住一个原则：取模型中部固定窗口的空间平均，不要取全流域。颗粒重排、入口效应和出口回流都会污染统计值，中部窗口是唯一能同时避开这三类干扰的位置。

本文还有配套的精品资源，点击获取