企业官网建设流程全解析

简介：一套基于MATLAB实现的格子Boltzmann方法两相流体动力学模拟系统，面向流体力学研究者、仿真工程师及数值计算方向学生，提供从模型构建到结果可视化的完整参考实现，解决气液两相流界面追踪与相变模拟中编程实现难度高的问题。压缩包共含2个文件，以m脚本和md文档为主：main.m覆盖格子Boltzmann方程离散化、碰撞流动迭代及两相界面追踪等核心流程，README.md说明模型参数设置、边界条件与运行方式，整体仅4KB，结构轻量，便于按需修改。目前已有82人学习，适合希望快速理解LBM两相流算法主干的MATLAB用户作为起点。通过阅读源码可掌握离散化处理、参数耦合和可视化输出的完整链条，减少从零搭建模型的时间成本；配套说明文档还能帮助排查边界条件与收敛性设置中的常见问题，为后续扩展到化工反应器、海洋工程或油水混合等场景打下基础。

1. 基于MATLAB的格子Boltzmann方法两相流体动力学模拟：先看现象，再谈原理

用MATLAB写格子Boltzmann方法的两相流体动力学模拟，第一反应往往是“这能跑得动吗”。实际当你需要快速验证液滴合并、界面演化或毛细现象时，这套系统反而是最直接的工具：不需要求解压力泊松方程，界面由伪势自动涌现，几十行代码就能观察到相分离。它为三类人准备：要复现论文数据的初学者、做参数扫描的工程师、以及想在开题前先看到动态现象的课题组学生。下面我们从算法骨架说起，一路讲到能复现的最小脚本和最容易翻车的几个参数。

2. 算法骨架：D2Q9、BGK碰撞与Shan-Chen伪势怎么组合

2.1 D2Q9离散速度与BGK松弛：两相模拟的物理起点

格子Boltzmann方法不直接跟踪分子，而是在每个格点上保存9个方向的粒子分布函数 f_q。宏观密度和速度由各方向分布函数求和得到，演化分两步：碰撞——分布函数向平衡态松弛；迁移——分布函数沿离散速度方向移动到相邻格点。两相模拟里的界面力、表面张力都依赖离散速度集的方向对称性，所以D2Q9是最常用起点，D2Q5只有正交方向，界面力的各向同性不足，会出现明显的网格痕迹。

D2Q9的速度配置和权重如下，记住权重之和等于1，声速 c_s² = 1/3，这在MATLAB实现里就是一组数组。

BGK碰撞的核心是松弛时间 tau，它通过 nu = c_s²(tau - 0.5) 决定动力学粘度。tau 必须大于0.5，否则粘度是负的，算不了几步就会发散。很多初学者的第一个NaN就是在这里踩的：物理粘度换算到格子单位后，回头算 tau 竟然小于0.5。

平衡态分布函数对低速流动是：

f_eq,q = w_q * rho * (1 + 3*(e·u) + 4.5*(e·u)^2 - 1.5*u·u)

这段公式在MATLAB里用数组乘法一次性算完，不要把格点写成双层for循环。

2.2 界面力从哪来：伪势模型的公式与局限

两相界面不是用跟踪网格或Level Set画出来的。Shan-Chen伪势模型在每个格点上定义一个伪势 psi(rho)，然后算格点与邻居的相互作用力：

F(x) = -G * psi(x) * Σ_q w_q * psi(x + e_q) * e_q

G是耦合强度，负值代表粒子间相互吸引。液相格点密度高，psi大，会吸引周围格点，于是高密度区自动聚成一团，低密度区成为汽相，界面就涌现出来了。最经典的伪势选择是：

psi = rho0 * (1 - exp(-rho / rho0))

rho0是参考密度，一般取1.0。这个形式的优点是数值上稳定，缺点是能做到的密度比有限，通常液相汽相密度比到5~20就接近极限。想要更高的密度比，需要换成基于状态方程的伪势，常见做法是用Carnahan-Starling状态方程反推psi。我的建议是先把指数形式的简单版本跑通，再考虑升级EOS版本。

2.3 最小可运行骨架：初始化、碰撞、迁移三步

这一节给出系统的最小骨架，真正的完整循环在下一章展开。先用数组把D2Q9配置写好，碰撞部分直接向量化，迁移用circshift。

% D2Q9 配置 w = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; cx = [0, 1, 0, -1, 0, 1, -1, -1, 1]; cy = [0, 0, 1, 0, -1, 1, 1, -1, -1]; cs2 = 1/3; % 假设 rho、ux、uy 已初始化，f 为 Ny×Nx×9 的数组 % 碰撞：BGK 松弛 for q = 1:9 cu = cx(q) * ux + cy(q) * uy; feq = w(q) * rho .* (1 + 3*cu + 4.5*cu.^2 - 1.5*(ux.^2 + uy.^2)); f(:,:,q) = f(:,:,q) - (f(:,:,q) - feq) / tau; end % 迁移：circshift 实现周期边界 for q = 1:9 f(:,:,q) = circshift(f(:,:,q), [-cy(q), cx(q)]); end

碰撞和迁移都只用9次循环，每次循环处理整个二维网格，这是MATLAB里效率最高的写法。注意circshift的第二个参数是[行偏移, 列偏移]，我把cy取负号，代表y正方向映射到行索引减小的方向，也就是视方向上移。如果你更习惯y轴向下增长，把这个符号整体对调即可，但迁移和后面力计算的符号必须保持一致，这是最容易出系统性错误的地方。

3. 从骨架到可跑系统：初始化、边界与主循环的落地细节

3.1 圆形液滴初始化：密度布点与平衡态赋值

初始化分两步：先给密度场布一个圆形液滴，再用局部平衡态给9个分布函数赋初值。很多入门代码直接从均匀场开始跑，让相分离自己发生，这样不仅收敛慢，而且可能长出多个小液滴而不是一个稳定的大液滴。我习惯直接放一个初始圆形液滴，让体系只做松弛。

% 200×100 格点，中心放半径 25 的液滴 Lx = 200; Ly = 100; rho0 = 1.0; radius = 25; [X, Y] = meshgrid((1:Lx) - Lx/2, (1:Ly) - Ly/2); rho = rho0 * ones(Ly, Lx); rho(X.^2 + Y.^2 <= radius^2) = 4.0 * rho0; ux = zeros(Ly, Lx); uy = zeros(Ly, Lx); f = zeros(Ly, Lx, 9); for q = 1:9 cu = cx(q) * ux + cy(q) * uy; f(:,:,q) = w(q) * rho .* (1 + 3*cu + 4.5*cu.^2 - 1.5*(ux.^2 + uy.^2)); end

初始速度设为零，所以平衡态公式里cu那一项就是0，但代码保留通用写法，后面加驱动或扰动时不用改初始化。液滴初始密度取4倍rho0，不是一个严格的稳定值，体系会在前几百步里自动调整到共存密度，这是正常的，不要因为初始密度变化就以为是发散。

meshgrid生成的X、Y坐标以液滴中心为原点，判断条件用X.^2 + Y.^2 <= radius^2，边界格点被包含进去，避免出现半个格点的锯齿。

3.2 边界条件怎么选：周期、反弹与压力边界的取舍

两相液滴在自由空间里扩散演化，最常见的选择是四周周期边界，circshift天然实现。但如果要做接触角研究，底部必须换成固体壁面，通常用半步反弹边界，把指向壁内的分布函数反向反射回去。

% 周期边界：迁移用 circshift，无需额外代码 % 上下壁面反弹示意：只写垂直两个方向，完整版需处理4个对角线分量 f(:, 1, 4) = f(:, 2, 2); % 底部壁面：向下分量换成来自上方邻居的向上分量 f(:,end, 2) = f(:,end-1, 4); % 顶部壁面对称处理

反弹边界的坑在于：只有垂直分量不够，对角线方向也要按镜像关系交换，否则壁面附近会出现速度泄漏。做壁面润湿时还需要在壁面格点放一个固定的伪势值，参与邻居求和，接触角大小就由这个壁面伪势决定。

压力边界在两相模拟里用得少，因为你需要同时给密度和速度约束，而两相体系的界面张力会使边界附近的密度场受干扰。如果只是跑液滴演化和Laplace校验，周期边界是完全够用的。

3.3 主循环与守恒监控：一个可复制的完整脚本

把前面的片段拼成一个完整主循环，这是我实际在用的模板：每一步按“宏观量计算 → 伪势力 → 碰撞 → 迁移 → 监控”的顺序执行。监控输出放在循环外，用mod控制频率，避免每步打印拖慢速度。

% 完整主循环：参数全部为格子单位 Lx = 200; Ly = 100; rho0 = 1.0; tau = 0.8; G = -5.0; radius = 25; nStep = 20000; nOut = 500; % D2Q9 配置 w = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; cx = [0, 1, 0, -1, 0, 1, -1, -1, 1]; cy = [0, 0, 1, 0, -1, 1, 1, -1, -1]; % 初始化 [X, Y] = meshgrid((1:Lx) - Lx/2, (1:Ly) - Ly/2); rho = rho0 * ones(Ly, Lx); rho(X.^2 + Y.^2 <= radius^2) = 4.0 * rho0; ux = zeros(Ly, Lx); uy = zeros(Ly, Lx); f = zeros(Ly, Lx, 9); for q = 1:9 cu = cx(q) * ux + cy(q) * uy; f(:,:,q) = w(q) * rho .* (1 + 3*cu + 4.5*cu.^2 - 1.5*(ux.^2 + uy.^2)); end sumRho0 = sum(rho(:)); for step = 1:nStep % 1. 宏观量 rho = sum(f, 3); ux = zeros(Ly, Lx); uy = zeros(Ly, Lx); for q = 1:9 ux = ux + cx(q) * f(:,:,q); uy = uy + cy(q) * f(:,:,q); end ux = ux ./ rho; uy = uy ./ rho; % 2. 伪势与作用力 psi = rho0 * (1 - exp(-rho / rho0)); Fx = zeros(Ly, Lx); Fy = zeros(Ly, Lx); for q = 1:9 psiN = circshift(psi, [-cy(q), cx(q)]); Fx = Fx - G * w(q) * psi .* psiN * cx(q); Fy = Fy - G * w(q) * psi .* psiN * cy(q); end % 3. 碰撞，作用力折入平衡态速度 uxEq = ux + tau .* Fx ./ rho; uyEq = uy + tau .* Fy ./ rho; for q = 1:9 cu = cx(q) * uxEq + cy(q) * uyEq; feq = w(q) * rho .* (1 + 3*cu + 4.5*cu.^2 - 1.5*(uxEq.^2 + uyEq.^2)); f(:,:,q) = f(:,:,q) - (f(:,:,q) - feq) / tau; end % 4. 迁移，周期边界 for q = 1:9 f(:,:,q) = circshift(f(:,:,q), [-cy(q), cx(q)]); end % 5. 监控 if mod(step, nOut) == 0 umax = max(max(sqrt(ux.^2 + uy.^2))); fprintf('step=%6d mass_err=%.3e u_max=%.3e\n', ... step, sum(rho(:)) / sumRho0 - 1, umax); end end

注意，宏观速度ux、uy是在碰撞之前从分布函数算出来的，而碰撞用的平衡态速度是uxEq、uyEq，这就是把伪势力折入动量方程的标准做法。伪势力用9邻域权重加权邻居的psi，和碰撞共用同一份w和cx、cy数组，符号也必须一致。

sum(rho(:))这种写法兼容旧版本MATLAB，不建议用sum(rho,'all')，后者要R2018b以上才支持，课题组里有人用旧版时很容易被这个函数卡住。

3.4 参数表与量纲说明

整套系统所有量都是格子单位：格距dx=1，时间步dt=1，松弛时间tau无量纲。物理粘度和速度都要用特征尺度自己换算，我建议把换算函数单独写成一个模块。

这套参数在200×100网格上跑2万步，通常5000步内界面就稳定了。u_max会缓慢下降到10^-3量级，这是伪速度，不是真空速度，不要用普通CFD的标准来判断收敛。

4. 调参指南：tau、G、密度比与界面厚度的关联

4.1 tau与粘度：先保住数值稳定再谈物理换算

tau直接决定数值稳定性。在BGK模型里，tau越接近0.5，粘度越低，但高频数值振荡越严重，界面附近的伪速度也越大。我的一般做法是先固定tau=0.8跑通，确认现象合理，再往低调。不要一开始就追求低粘度，那会把“参数问题”和“算法问题”混在一起，排错难度翻倍。

从物理粘度反推tau的公式是：

tau = nu_lattice / cs2 + 0.5 其中 nu_lattice = nu_physical * dt / (dx^2)

我见过最典型的翻车案例是把水的10^-6 m²/s直接当nu_lattice代入，算出tau约等于0.5，跑几十步就NaN。因为格子单位下的dx、dt必须先通过特征尺度确定，物理粘度除以dx²再乘dt，量级完全不一样。

4.2 G与伪势：界面厚度、密度比和稳定性怎么取平衡

G的绝对值决定了“吸引力”的强弱。G太弱，界面带厚，两三格内密度从液相过渡到汽相，表面张力效果不明显；G太强，力项导致密度振荡，液滴边缘出现负密度。做参数扫描时，我习惯先固定tau和rho0，从G=-2开始，每步加1，观察界面密度剖面和u_max。

% 参数扫描模板：G从-2扫到-6 GList = -2:-1:-6; for k = 1:numel(GList) G = GList(k); % 重新初始化 rho 和 f % 运行到稳态后记录密度比、界面厚度和最大伪速度 rhoLiq = mean(rho(rho > 2.0)); % 液相平均密度 rhoVap = mean(rho(rho < 0.8)); % 汽相平均密度 fprintf('G=%g rho_ratio=%.2f\n', G, rhoLiq/rhoVap); end

液相和汽相的阈值不是固定的，我这里是针对rho0=1、初始液滴密度4倍的设置。换rho0后要重新看密度分布直方图确定阈值。密度比的限制不是bug，是指数型伪势模型本身的特性，想要更高的密度比就去换EOS型伪势，那是另一个量级的调试工作量。

4.3 监控指标：用界面带速度和总质量判断收敛

两相模拟的收敛判据，看全局动量没意义，因为体系最终趋于静止。我习惯两个指标：界面带平均伪速度和总质量相对误差。

% 界面带：密度介于液相和汽相之间的格点 mask = abs(rho - (rhoLiq + rhoVap)/2) < (rhoLiq - rhoVap)*0.15; uSpur = sqrt(mean(ux(mask).^2 + uy(mask).^2)); if uSpur < 1e-6 fprintf('converged at step %d\n', step); break; end

界面带的选取决定了伪速度的量级，阈值0.15对应15%的界面宽度，这是我试下来比较稳的设置。uSpur能压到10^-6以下，通常说明已经平衡。如果uSpur停在10^-4不再下降，不要继续加时间步，先回去调G和界面厚度。

4.4 无量纲参数对照表

给结果配图或写入报告时，别人问的第一句话往往是“这个仿真对应多少雷诺数”。我建议把关键无量纲组列成表格，报告里直接复用。

换算表面张力sigma用Laplace定律标定，这一步单独在后面讲。调参阶段不要纠结sigma准不准，先用一个固定值把界面和流动趋势跑出来。

5. 避坑：发散、界面断裂与质量漂移的排查记录

5.1 跑几十步就NaN：先查tau再查G

现象：程序在100~200步内输出NaN，时间早于界面稳定。

原因：九成是tau小于等于0.5，粘度失去物理意义；剩余一成是G绝对值太大，伪势力造成的密度振荡直接击穿rho。还有一个隐蔽来源：宏观速度除rho时，rho先出现负值，之后所有除法都崩掉。

解决：先用一行代码检查参数合法性。if tau <= 0.5 || G >= 0直接报错。再把G降到-3以下重新跑。如果仍然NaN，在每一步末尾打印min(rho(:))，看到负数出现在液滴界面附近，说明是力项太强，而不是tau问题。我自己的调试顺序永远是：先查tau，再查G，最后查符号。

5.2 液滴界面偏移或断裂：迁移方向与力方向不一致

现象：初始圆形液滴演化后变成椭圆，或界面沿某个轴出现锯齿，整体形状不对称。

原因：迁移用的circshift符号和伪势力里的邻居求和符号不一致。力计算里psiN取邻居psi，方向用cx、cy；迁移里分布函数移动方向也用cx、cy，如果两个地方一个取负一个不取负，物理上等于力向右推、粒子向左跑。

解决：把符号统一成一套。我最开始写的时候也翻过车，排查方法是在界面稳定后输出速度场矢量图，如果速度场整体朝一个方向飘，说明力方向反了；如果速度场呈对称的“涡对”，则符号一致但G方向有问题。

% 错误写法：迁移的cy符号与伪势力不一致 % f(:,:,q) = circshift(f(:,:,q), [cy(q), cx(q)]); % 正确写法 % f(:,:,q) = circshift(f(:,:,q), [-cy(q), cx(q)]);

5.3 总质量缓慢漂移：边界不对称或人为截断

现象：mass_err在几百步后从10^-14逐渐涨到10^-6甚至更高。

原因：周期边界下circshift本身严格守恒，漂移几乎都来自“人为干预”。常见有三处：一是在碰撞后对rho做了钳制，比如rho = max(rho, 0)；二是边界用了反弹但只处理了垂直方向，对角线分布函数在壁面泄漏；三是用单精度后监控方式不对，把累积误差当成质量损失。

解决：不要在碰撞循环里修改rho，负密度出现在参数失控之后，而不是解决手段。反弹边界要写完整版，四个对角线方向按镜像交换。监控误差用相对值，单精度下10^-6以内的漂移属于正常，不要追求机器零。

5.4 伪速度降不下去：先确认是不是模型本身的限制

现象：界面带最大速度u_max停在10^-4~10^-3量级，怎么调tau都没用。这个现象在Shan-Chen模型里是模型特性，不是数值bug。

原因：离散格点的各项异性、界面厚度只有几个格点、密度梯度在界面处非常大，三者共同产生了所谓的spurious currents。界面厚度越小、密度比越高，伪速度越大。

解决：先接受它，报告里写清楚格子单位下的u_max即可。如果应用场景不能容忍伪速度，有三个方向：增大界面厚度，把G绝对值调小；改用精确差分法（EDM）替代速度修正法折入力；换多层分布函数或多松弛时间(MRT)模型。其中EDM是我见过最有效也相对容易实现的升级，改动只在力折入平衡态那一步。

6. 验证这套系统的两板斧：Laplace定律校验与性能优化

6.1 用Laplace定律校准表面张力

跑通液滴演化后，第一件该做的事是验证表面张力：对一组不同半径的液滴，分别跑到稳态，记录液滴内外压力差，做dp对1/R的线性拟合。线性关系成立，说明界面力计算是各向同性的；斜率就是表面张力。

% 每个半径跑出稳态后，提取内外压差 maskL = rho > 2.5; % 液滴区 maskV = rho < 1.5; % 汽相区 p = 1/3 * rho; % 相对压力，用于线性关系验证 dp(k) = mean(p(maskL)) - mean(p(maskV)); % 最后对 dp 与 1/radius 做线性拟合，斜率即 sigma

直接拿1/3*rho当压力只适合做相对比较。伪势模型里压力项还包含相互作用贡献，严格校准要用带伪势修正的状态方程。我的习惯是先用这个简化值确认dp与1/R线性，再看具体的sigma是否落在预期量纲范围；若要精确值，再补上EOS修正项。

6.2 接触角与壁面润湿：一个半小时内能出结果

把周期边界换成底部反弹壁面，在壁面格点放固定伪势，就能测接触角。稳态后提取液滴底宽w和顶点高度h，按球冠假设反推接触角：

% 从二值化轮廓提取 w 和 h 后 theta = 2 * atan(2 * h / w); % 弧度，乘以180/pi转角度

测量时注意两点：一是取平衡态帧，不要取演化中途的轮廓；二是接触角对壁面伪势的初值很敏感，先扫描一组壁面伪势值，把接触角-伪势曲线做出来，再选择目标角度对应的参数。这是两相模拟里最容易出漂亮图表、也最容易被审稿人追问的验证手段。

6.3 再快一点的优化习惯

MATLAB跑LBM，最大的瓶颈不是exp函数，而是9次circshift的数组拷贝。三个实用习惯：整个程序用single精度存储f，速度提升明显，监控值不受影响；参数扫描用parfor，每个G值独立跑；对大网格，先把索引表算好，迁移时用索引代替每次circshift。

% 把迁移写成索引表，减少重复计算 idxX = circshift((1:Lx)', 1); idxY = circshift((1:Ly)', -1); % 迁移时直接用 f(:,:,q) = f(idxY, idxX, q)

我最初做这个系统时，把大量时间花在调G和追界面断裂上，后来才明白这属于方法本身固有的参数窗口，不是代码写错。真正值得花时间的，是先把D2Q9骨架和监控指标固定下来，再谈模型升级。希望这些排查顺序和参数边界能帮到你。

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