多相流听起来是热能与动力工程、化工机械的专属课题但只要你在烧烤摊前盯着啤酒里上升的气泡出过神就已经接触过它的本质。多相流程序这个名词如果放在MATLAB语境下就是让电脑替你“看”气泡怎么上升、颗粒怎么翻滚、油和水怎么在管道里你追我赶。作为一个从写脚本小白一路踩坑走过来的MATLAB用户我想把这些经验整理成一篇真正能照着做的文章而不是教科书式地贴公式。这篇文章适合刚接触多相流模拟的研究生、想做方案预研的工程师以及想搞懂数值求解底层逻辑的爱好者。你会知道MATLAB为什么适合干这件事、数学模型怎么落地成代码、以及最常见的数值爆炸怎么避免。1. 多相流问题为什么值得用MATLAB来处理1.1 多相流建模的核心矛盾尺度与精度多相流的本质是两种及以上相态气、液、固在同一个空间里存在并且相互之间不断传质、传热、交换动量。范围可以从化工反应器里的气液鼓泡塔到石油输送管道里的油气水三相混合物再到流化床里数百万颗粒的翻滚。不同场景下你关注的物理量完全不一样有时候只要管道总压降有时候必须知道气泡的尺寸分布有时候甚至要捕捉液滴和界面之间那层薄膜。这就带来第一个核心矛盾尺度与精度不可兼得。如果你把两相看成一种均匀混合物计算极其省事但气泡在哪里、滑移速度多大完全说不清如果你把每一个气泡、每一颗颗粒都单独建模精度上去了计算量却呈指数级膨胀。所以写多相流程序的第一步根本不是打开MATLAB而是想清楚你的问题需要什么尺度下的什么精度。我见过太多人一上来就想复现文献里的三维气泡群结果困在网格和收敛里好几个月其实他手里的课题用一维漂移流模型就能给出足够好的结果。1.2 MATLAB相对专业CFD软件的优势与边界过去我在课题里同时用过FLUENT和OpenFOAM也用过MATLAB自己写程序。说实话专业CFD软件在网格划分、湍流模型、大规模并行上确实强但多数的多相流研究阶段是探索性的模型要反复改参数要来回试可能还要跟实验数据一遍遍比对。这时候MATLAB的优势就很明显了。首先是矩阵运算天然适配离散后的代数方程组你不需要去点一堆GUI按钮所有求解逻辑都摊在眼前出错能迅速定位。其次是可视化非常方便pcolor、quiver、contour几个命令一组合就能把相场和速度场画成动画做学术汇报、写论文插图都很顺手。但MATLAB也有很清晰的边界。它不适合做几百万网格的三维工业级模拟纯解释型循环在大规模计算上会慢到让你怀疑人生它也没有内置的全功能多相流CFD求解器至少没有一个能开箱即用的“双欧拉VOF群体平衡”全家桶。所以我的定位一直是MATLAB适合做机理研究、方案预研、教学演示以及小规模的二维或轴对称问题。如果你最终需要极其精细的工程级结果可以用MATLAB搭模型验证算法再移植到专业CFD工具里。这不是绕路反而是最稳妥的路线。2. 从物理模型到可执行的MATLAB程序2.1 三种常用多相流数学模型的适用场景多相流数学模型多到能让人眼花缭乱但核心逃不开下面三种框架。模型基本思想适用场景编程难度均相流把多相当成一种均匀混合物速度和温度取平均低含液率气体的管道压降估算、热力循环简化低漂移流在均相基础上加入相间滑移用漂移速度修正竖直管中的气液两相流、泡状流、池沸腾中两流体模型每一相分别建立守恒方程通过相间作用力耦合流型复杂、需要相分布和局部滑移的精细问题高我个人建议初学者先把均相流和漂移流吃透。它们不需要处理复杂界面MATLAB里几十行脚本就能得到很有物理意义的结果。比如竖直管内气液两相流你用漂移流模型算出截面含气率再用摩擦压降关联式求出压降跟实验数据对比往往趋势非常漂亮。这个过程会帮你建立“程序不过是对物理方程的离散翻译”这个直觉。2.2 对偏微分方程组做时空离散有限差分与有限体积多相流控制方程大多是偏微分方程组特别是连续性方程和动量方程。MATLAB不会直接帮你解偏微分方程真正要你做的是把连续的偏微分方程在网格上离散成代数方程。两种主流思路是有限差分和有限体积。有限差分实现简单直接对导数用差分公式近似有限体积则把每个网格当作控制容积积分守恒性更好处理间断比如相界面时更不容易出幺蛾子。我写最基础的平流方程练手时会这样实现迎风差分。假设要计算一个一维方波在速度场中的输运代码如下L 1; Nx 200; dx L / Nx; x (0:Nx-1) * dx; v zeros(1, Nx); v(1:30) 1; % 初始方波 u 0.5; % 对流速度 dt 0.001; nsteps 500; for n 1:nsteps vn v; for i 2:Nx-1 if u 0 v(i) vn(i) - u * dt / dx * (vn(i) - vn(i-1)); else v(i) vn(i) - u * dt / dx * (vn(i1) - vn(i)); end end % 简单边界条件 v(1) v(2); v(Nx) v(Nx-1); end % plot(x, v); xlim([0 L]); title(迎风差分后的方波);注意这里用了一个关键选择迎风格式。为什么不用中心差分因为中心差分对对流项天然会产生振荡而迎风格式依据流动方向取上游值虽然有一阶耗散界面会被抹平一点但至少数值上是稳定的怪物不会爆炸。初学者往往没意识到格式的选择不是玄学而是和物理过程一一对应的。2.3 压力-速度耦合与CFL条件的通俗解释一旦动量方程里出现压力梯度程序就会面临压力-速度耦合问题。你可以把它想成一群人在排队进地铁速度的改变取决于后面气流推力压力梯度但推力大小又反过来取决于人群密度速度散度。最简单的处理办法是人工可压缩法或投影法。投影法的核心是先求一个临时速度再用压力修正让它满足连续性方程相当于先跑一步再拉回来。这个流程在MATLAB里其实不复杂只是需要在每个时间步内做一次或几次迭代。另一个绕不开的稳定性判据是CFL条件。信息传播的速度乘以时间步长不能超过一个网格尺寸说得直白点你骑着共享单车车时速20公里每秒钟至少要在一个路灯之间刷新一次如果总时间步太大就会一路跳过好几个路灯完全丢失中间信息。公式写作dt dx / max(|u|)如果dx 0.01 m最大速度是1 m/s那dt至少小于0.01 s。我习惯取0.5倍的临界值也就是0.005 s。这个习惯看起来保守但能省下一大堆排查发散的时间。3. 亲手写一个气泡上升模拟程序二维VOF简化版3.1 程序框架与网格初始化理解了离散概念就可以动手做一件微妙的事模拟一个气泡在液体中上升。我选择二维矩形计算域80格宽、160格高初始在底部偏上位置放置一个半径15格的圆形气泡。为了让问题简单先忽略表面张力、只考虑浮力和粘性力并使用VOF流体体积函数方法表示两相用一个标量alpha代表液体体积分数alpha等于1就是纯液体等于0就是纯气体在0到1之间就是界面区域。初始化代码如下思路就是先铺网格再根据几何位置给alpha赋值nx 80; ny 160; dx 0.001; dy 0.001; % 无量纲网格尺寸 x (0:nx-1) * dx; y (0:ny-1) * dy; alpha zeros(ny, nx); % 1为液体0为气体 center_x 0.04; center_y 0.04; R 0.015; rho_l 1000; rho_g 1.2; mu_l 0.001; mu_g 1.8e-5; rho zeros(ny, nx); mu zeros(ny, nx); for i 1:nx for j 1:ny dist sqrt((x(i)-center_x)^2 (y(j)-center_y)^2); if dist R alpha(j,i) 0; else alpha(j,i) 1; end rho(j,i) alpha(j,i) * rho_l (1-alpha(j,i)) * rho_g; mu(j,i) alpha(j,i) * mu_l (1-alpha(j,i)) * mu_g; end end imagesc(x, y, alpha); axis xy; colormap([1 1 1; 0 0 0]); colorbar;这段代码里另外生成了rho和mu场。很多初学者一开始只追踪alpha不考虑密度和粘度的空间变化结果浮力/粘性力算出来全是错的。记住多相流里的物性场必须跟随alpha实时更新。3.2 相场输运与界面重构的核心步骤气泡上升的核心方程是alpha的输运方程∂α/∂t ∇·(α U) 0我们必须在每个时间步求解这个方程然后重新计算密度、粘度场再去解流场。这一步的数值格式决定了界面到底保持得多锋锐。只用简单的中心差分界面很快会糊成一片灰色过渡带气泡看起来像一团雾。我一般会引入一个“人工界面压缩”项在界面附近加入一个和法向相关的反扩散项用来抵消数值抹平。实际上就是用迎风格式加一个压缩通量。简化后的核心循环长这样% 假设已经计算得到速度场 u_xi, u_yi % 求解 alpha 输运 alpha_new alpha; for i 2:nx-1 for j 2:ny-1 flux_x (u_xi(j,i) 0) * alpha(j,i-1) (u_xi(j,i) 0) * alpha(j,i); flux_y (u_yi(j,i) 0) * alpha(j-1,i) (u_yi(j,i) 0) * alpha(j,i); alpha_new(j,i) alpha(j,i) - dt/dx*(u_xi(j,i)*alphax) ... - dt/dy*(u_yi(j,i)*alphay); end end实际写出正确通量需要小心处理迎风这里不展开完整代码。更关键的是每个时间步要把alpha截断到[0,1]之间否则质量守恒会崩。你可以把它想象成相机曝光过度稍微超出一点没关系超出太多画面就全白了。程序里我用一句alpha_new max(0, min(1, alpha_new));来兜底。3.3 可视化与动画输出MATLAB最让我喜欢的就是可视化随手就来。展示alpha场可以用imagesc或pcolor叠加速度场用quiver。如果你要生成论文里那种动态效果最好录制视频。基本套路是v VideoWriter(bubble_rise.avi); open(v); for it 1:nsteps % ... 计算更新 ... if mod(it, 5) 0 figure(1); imagesc(x, y, alpha); axis xy; title([t , num2str(it*dt)]); hold on; % 画流线或速度矢量 ux_plot u_xi(1:4:end, 1:4:end); uy_plot u_yi(1:4:end, 1:4:end); quiver(x(1:4:end)/dx, y(1:4:end)/dy, ux_plot, uy_plot, k); hold off; axis equal; xlim([0 nx*dx]); ylim([0 ny*dy]); drawnow; writeVideo(v, getframe(gcf)); end end close(v);有一个小坑pcolor默认的网格显示会有白线空隙最好加上shading interp。同时注意坐标轴纵横比气泡是圆的别因为x和y的物理尺寸不一致给拉成椭圆。3.4 结果解读与验证看什么指标程序跑通之后最忌讳的事情就是盯着彩色云图喊“哇好漂亮”然后直接写论文。你必须量化验证。对于气泡上升我常看三个指标气泡重心位置随时间的变化用来计算上升速度、气泡等效直径的演化、总液体体积分数的守恒理论上误差应该低于几个百分点。气泡重心位置很简单每个时刻对气泡区域alpha0.5的网格计算面积加权平均的y坐标然后求差除以时间步。算出来的上升速度如果接近理论终端速度量级说明基本物理对。我的经验是如果你的二维气泡终端速度比相同直径的三维气泡大一些也别慌因为二维和三维的绕流阻力本身就有差异这种由维度引起的差异可以用文献数据定性对比不必追求完全吻合。4. 颗粒-流体两相流一个简单的离散元DEM耦合案例4.1 为什么颗粒流要考虑接触力与拖曳力除了气泡这类界面问题另一种常见多相流是颗粒伴随流体运动典型代表就是流化床和气力输送。这时候一个颗粒受到的力包括重力、浮力、流体拖曳力和与其他颗粒的接触力。拖曳力模型有太多我常用Wen-Yu或者Ergun模型核心是拖曳系数β与空隙率、颗粒雷诺数相关。接触力则简化为弹簧阻尼模型法向力正比于颗粒重叠量阻尼项耗散碰撞能量。这些力加起来就构成颗粒的牛顿第二定律。为了演示方便可以先做“流体静止、颗粒沉降”的简化案例一个颗粒在液体中自由沉降计算拖曳力、浮力、重力更新位置。加上多个颗粒后才需要处理碰撞和空隙率统计。这个渐进过程能避免一上来就被复杂的耦合搞晕。4.2 MATLAB实现的时间推进与邻居搜索DEM程序的时间推进思路跟求解流体类似每个时间步计算颗粒受到的合力然后更新速度和位置。但碰撞力要求知道颗粒之间是否接触所以需要邻居搜索。小规模颗粒几百个以内全对搜索也能跑Np 200; rad 0.0005; % 半径 pos rand(Np,2) * 0.02; vel zeros(Np,2); dt 1e-5; k 1000; c 10; m 0.001; g 9.8; for t 1:1000 force zeros(Np,2); force(:,2) -m*g; % 重力 % 与流体拖曳力等略这里只演示接触 for i 1:Np for j i1:Np dx pos(i,1)-pos(j,1); dy pos(i,2)-pos(j,2); d sqrt(dx^2dy^2); if d 2*rad delta 2*rad - d; nx dx/d; ny dy/d; fi - (k*delta c*dot(vel(i,:)-vel(j,:), [nx ny])) * [nx ny]; force(i,:) force(i,:) fi; force(j,:) force(j,:) - fi; end end end vel vel (force/m) * dt; pos pos vel * dt; pos(posrad) rad; % 简单边界 pos(pos0.02-rad) 0.02-rad; end这个示例省略了流体耦合但已经能看出DEM的最小闭环。真正与CFD耦合时每步还要把颗粒位置映射到流体网格计算空隙率再把拖曳力反馈给流体动量方程。你会发现“耦合”的本质就是在流体场和颗粒场之间反复传递源项。4.3 如何稳定追踪颗粒并统计空隙率写DEM最痛苦的就是时间步长。颗粒碰撞时间尺度通常远小于流体平流的时间尺度所以我一般取dt小于颗粒接触时间为一个数量级以上。否则颗粒会像子弹一样互相穿透场面看起来不像是物理模拟而是物理爆炸。空隙率统计常用“中心点法”或“体积分数法”。体积分数法是算每个网格内颗粒体积占比然后用1减去它就得到空隙率ε。一个容易踩的坑是当网格尺寸小于颗粒直径时一个颗粒可能跨好几个网格统计需要按体积加权。如果你只是粗略看流化床的整体空隙率中心点法其实够用了。我建议在小规模练习里先画颗粒位置散点图再把空隙率云图画成背景这样能直观看到流体通道结构。5. 数据后处理与经验避坑清单5.1 用MATLAB对模拟结果做统计和可视化多相流程序跑出来的原始数据往往是一堆矩阵和数组不后处理根本没法用。我常用的策略是写一个后处理脚本把不同时刻的结果打包统计。比如气泡尺寸分布可以用histogram实现颗粒轨迹用scatter叠加颜色映射速度沿管道径向统计含气率分布则可以用三维数组的切片操作配合mean。还有一类重要后处理是与实验数据对比把模拟得到的压降、含气率等整理成表格用lsqnonlin拟合经验公式里的参数。这个环节常常能反馈出模型选错或边界条件不对的问题。5.2 多相流程序常见错误和排查方法速查表在实际运行中大部分时间不是花在写新功能而是花在排查程序为什么发散、界面为什么糊、颗粒为什么飞。下面这张表是我踩坑后的浓缩版。症状可能原因解决思路NaN/Inf时间步过大、密度比异常、边界条件不稳定降低dt检查物性场强制截断负体积分数迎风格式不满足、步长过大加通量限制器减小dt界面过于模糊数值扩散严重使用高阶格式或人工界面压缩项颗粒重叠或穿透接触时间步不够小将dt降到碰撞特征时间的1/10以下压力迭代发散初始速度场散度不为零初始化满足连续性使用投影法内存不足网格过密、变量存储过多用稀疏矩阵减少中间变量考虑GPU这里我想强调第一类问题NaN。有一次我模拟高密度比气液两相液体密度是气体的一千倍结果刚算几十步就全盘NaN。排查了半天发现是初始化时把气相网格的密度设成了0导致动量方程里出现无穷大。后来我在所有物性初始化后加了一行检查if any(rho(:)0) error(density positive); end省了无数后续麻烦。5.3 提高效率与稳定性的4个微调技巧最后聊聊让程序跑得更快、更稳的经验。第一个技巧是自适应时间步。不要傻傻地固定dt每步根据当前最大速度重新计算cfl 0.5; dt cfl * min(dx, dy) / max(max(abs(u_xi)), max(abs(u_yi)));这能把速度慢时的大量冗余步数省掉也在速度变快时自动保护稳定性。第二个技巧是通量限制器。迎风格式虽然稳但太“糊”格式阶数高点又容易振荡。折中方案是使用带限制器的线性格式比如minmod或superbee。MATLAB实现并不复杂很多材料里都能找到现成函数。第三个技巧是隐式处理粘性项。粘性项的显式稳定性条件里时间步和网格间距的平方成正比。一旦网格加密dt会小到离谱。把粘性项改成隐式或者半隐式就能让时间步只受对流限制计算量立刻下降一个档次。第四个技巧是并行扫描。你需要扫多组参数时可以用parfor同时跑不同工况。我上次模拟不同气泡直径的终端速度用parfor一次性把20组工况算完省下来一整个下午。6. 从单机脚本到系统级多相流仿真6.1 用Simulink搭建气液两相流系统动态模型有时你不关心气泡具体的运动轨迹只想研究整个系统动态比如分离器液位控制、管道气液两相流压力波动。这种时候逐格求解CFD不经济应该做集总参数模型把质量守恒和动量守恒简化成一组常微分方程。Simulink 里搭建这种系统很顺手用积分模块表示液位/压力状态用函数模块计算相间质量交换和动量源项然后拉一个PID控制器形成闭环。我见过很多做控制策略的研究者先用简化的Simulink两相流模型验证算法再拿到专业软件里跑三维细节效率非常高。6.2 与实验数据联动参数标定与模型校准多相流模型里经验关联式中的参数往往来自特定实验范围换一套工况就会失准。这时可以用MATLAB优化工具箱做参数标定。基本流程是把实验数据导入为表格写一个目标函数计算模型预测值与实验值的残差然后调用lsqnonlin或patternsearch寻找最优参数。重点是要把参数取值范围限制在物理合理区间否则优化器会把某个系数调到负几百纯粹变成数学游戏。我还会刻意保留一组“验证工况”不参与拟合用来检测模型有没有过拟合。6.3 后续扩展方向当你把一个二维气泡上升程序或者DEM与CFD耦合程序跑通后接下来的扩展方向其实很多。你可以加表面张力项用连续表面力模型替代目前的忽略处理也可以把二维数组扩展到三维并用gpuArray做加速更可以积累一批模拟结果用MATLAB的机器学习工具箱训练一个快速预测模型以后毫秒级输出终端速度或压降做工程筛选时非常实用。我个人体会是不要一上来就追三维和GPU先在二维把数值格式和物理过程吃透后续扩展会顺畅得多。最后再分享一个小技巧我电脑里始终存着一个“最简多相流模板”只包含网格生成、一个平流方程的更新、一张结果动画、以及几行防发散检查。每次写新项目都在这个模板上不断增加内容。多相流程序的奇妙之处就在于它把复杂的物理世界压缩成矩阵和循环又在屏幕上以瞬变的图像释放出来。当你第一次看到自己写的气泡从底部摇摇晃晃升起来那种成就感比看别人视频里的漂亮结果强一百倍。
