MATLAB多相流仿真从建模到代码：颗粒沉降与漂移通量模型实战

我是被这个问题找上门的有人要在MATLAB里写多相流程序模拟气泡在液体中的运动还要能输出动画和曲线。他搜了一圈资料要么是商业CFD软件的界面截图要么是几十年没人维护的Fortran老代码看完更懵了。这个场景我太熟悉了。MATLAB多相流程序并不是某个现成的工具箱而是一整套建模思想和数值方法的组合MATLAB只是把你脑子里那套物理模型快速变成可运行代码的试验场。这篇文章把我实际写多相流程序时的路线图整理出来从建模方法选型入手给了两个可以直接跑的案例一个是颗粒沉降的拉格朗日粒子追踪一个是垂直管内气液两相流的漂移通量模型。代码不长物理背景我会讲透参数也会逐个说明为什么这么取最后附上我踩过的那些坑。不管你是在做课程设计、研究生课题还是刚转行做多相流仿真照着这个思路走至少不会再对着空白编辑器发呆了。1. 多相流建模分类与MATLAB的适配边界1.1 多相流动到底在算什么多相流顾名思义是两种或多种相态同时流动的系统。日常生活里常见的气泡在水中上升、油水在管道中分层、流化床里颗粒被气流吹起来都是多相流。工业上更不陌生换热器里的沸腾、反应器里的气液搅拌、油气集输管线中的段塞流全属于这个领域。多相流难在什么地方因为不同相之间有界面界面会变形、移动、破碎、合并。单相流动只需要求解一个速度场和压力场多相流则还要追踪相的分布描述相间的相对滑移处理表面张力、曳力、相间传热传质。你没法用一套固定的网格来描述一个不断变化的界面这就是它比单相流复杂得多的根源。但如果从程序编写的角度想问题事情就简单一些。多相流程序要回答的核心问题其实只有三个各相在空间里怎么分布、各自以什么速度运动、相与相之间如何相互作用。不同的建模方法本质上是对这三个问题的不同回答方式。搞清楚了这一点再去看MATLAB里的代码实现脉络就清楚了。1.2 三大建模框架欧拉-欧拉、欧拉-拉格朗日与界面追踪我给大家一个非常不严谨但很好用的类比。假设你要研究一个城市广场上的人流欧拉-欧拉方法是把行人和空气都抽象成相互穿透的连续流体不看单个的人只看每一块区域里行人浓度和空气浓度随时间怎么变。这种方法适合相含率高的系统比如流化床、气泡塔、喷雾塔两相都是连续介质各自有速度和体积分数通过相间曳力耦合。工程上大量用但代价是丢失了界面细节你不可能看出某个气泡长什么样。欧拉-拉格朗日方法是连续相仍然用欧拉框架求解但离散的颗粒、液滴、气泡则被当成一个个质点追踪。这相当于你只关心广场上每个具体的人走哪条路不关心每一个人脚底下的空气流场细节。它很适合稀相系统颗粒体积分数一般小于10%比如喷雾干燥、气力输送、颗粒沉降。好处是能直接看到颗粒轨迹坏处是颗粒多了之后计算量爆炸。界面追踪方法比如VOF流体体积法和Level Set水平集法则是直接去捕捉相界面的位置。它像你在广场上空架了个摄像头一点一点盯住人流的轮廓线在哪里。这个方法适合模拟气泡上升、液滴撞击、射流破碎这类界面形态剧烈的流动。VOF在商用软件里最常用但需要求解界面输运方程做表面张力处理对网格和时间步很挑剔。把这三种框架落到MATLAB里各有各的实现形态。欧拉-拉格朗日在MATLAB里实现最顺手本身就是一堆粒子位置和速度的更新循环。欧拉-欧拉则大量化简为偏微分方程组再用有限差分或pdepe去解。界面追踪在MATLAB里也能写教学版只是要接受它跑不过商业软件的事实。建模框架离散相处理方式适合场景MATLAB实现难度典型代表欧拉-欧拉各相均为连续介质高含率工业设备中等需要解PDE双流体模型欧拉-拉格朗日粒子逐点追踪稀相颗粒/液滴低循环即可DPM粒子追踪界面追踪直接捕捉界面气泡/液滴/射流较高需输运方程VOF、Level Set1.3 MATLAB能干什么不能干什么先泼一盆冷水。MATLAB不是OpenFOAM也不是Fluent它没有成熟的非结构网格划分器没有并行大涡模拟能力也没有一整套湍流壁面函数库。你要拿MATLAB去模拟一个工业反应器内部的三维气液两相湍流那是拿错了工具。但MATLAB在多相流领域依然有它不可替代的位置。第一是教学和课程设计很多流体力学和数值方法课程需要学生亲手写一个简化求解器MATLAB语法友好调试方便比C门槛低太多。第二是科研验证当你提出一个模型修正式想快速看看趋势是否合理写个MATLAB脚本半天就能出结果。第三是算法原型在把代码移植到高性能平台之前先用MATLAB验证数值格式的正确性这是很多研究者实际在走的路。我的建议是如果目标只是看到一张漂亮的气泡上升云图请直接学OpenFOAM或者Fluent如果你的目标是搞清楚多相流模型内部每一步在算什么MATLAB是极好的地方。它逼着你把守恒方程、离散格式、迭代算法全部用代码写出来写一遍你就再也忘不掉了。2. 方案选型判断该用哪种多相流模型2.1 三个问题决定你的建模路线第一个问题你关注的到底是界面还是粒子想模拟液滴在气流中的变形和破碎你需要界面信息那就走VOF或Level Set路线。只想统计一群颗粒在管道里的浓度分布和停留时间那欧拉-拉格朗日就够了。第二个问题离散相浓度高不高浓度低到颗粒之间基本不碰撞拉格朗日单颗粒追踪就合理。浓度高了颗粒间相互作用和反向影响连续相都不能忽略这时候要么回到欧拉-欧拉双流体模型要么用更复杂的碰撞模型。第三个问题你是要一维趋势还是要二维/三维的细节很多时候工程上只需要一条压降曲线、一个空隙率变化趋势一维漂移通量模型几十行MATLAB代码就解决了根本不需要动用CFD求解器。先问自己这个问题能帮你省掉大量无意义的建模工作。我见过太多人一上来就说我要用VOF问为什么回答是因为别人都用VOF。但VOF在MATLAB里写完整教学版至少得几百行无数边界条件和表面张力细节要处理很多人直接就被磨掉耐心了。方案选型的本质是做减法在满足物理需求的前提下选最简单可靠的模型。2.2 拆解多相流问题的步骤我的习惯是拿到问题先做五个步骤不被具体物理场景带着跑第一步确定连续相和分散相。例如气泡在水中上升连续相是水分散相是气泡颗粒在水中沉降连续相是水分散相是颗粒。第二步估算分散相的尺寸、浓度和最重要的无量纲数。雷诺数Re判断流动是层流还是湍流斯托克斯数St判断颗粒跟随性韦伯数We判断界面是否容易破碎。这些无量纲数决定了你该用哪种阻力模型以及表面张力重不重要。第三步列出主导力。颗粒沉降里主要是重力、浮力和曳力气泡上升还有表面张力管道流动还要加摩擦压降和重力压降。不是所有力都重要先把主导项列出来。第四步选择建模框架。根据颗粒浓度和界面重要性用前面说的三个问题去判断。第五步寻找一个可以验证的极限情况。比如颗粒沉降的终端速度可以解析估算管道流动的空隙率可以用漂移通量公式直接算。有了这个标尺你代码写错了也能立刻发现。无量纲数的计算一定不要偷懒。我读模型代码第一件事就是找Re算没算怎么用的。很多人写粒子追踪阻力系数直接照抄斯托克斯公式完全不检查Re范围结果模拟出来终端速度偏差好几倍还不知道错在哪。2.3 MATLAB代码组织的三个建议第一把物理参数集中放在文件头部用清晰的变量名。不要在一个200行的脚本里到处散落魔法数字。我做过的所有靠谱的MATLAB仿真程序都遵守参数区、初始化区、主循环区、后处理区的四段式结构。第二预分配所有数组多用矩阵运算。MATLAB的循环性能差但如果你写的是粒子追踪每个粒子更新位置本质上就是一个矩阵操作可以一次性更新几千个粒子。后面颗粒沉降案例里我只写一个粒子但稍微改一两个下标就能扩展成粒子群。第三后处理早点做。不要等程序完全跑对了再画图。我会先画一条速度曲线看趋势对不对再继续往下调参数。用drawnow limitrate控制动画帧率既能看实时过程又不会被卡死。3. 案例一颗粒沉降——用MATLAB实现对颗粒轨迹的追踪3.1 物理模型与假设颗粒沉降是多相流里最经典的入门案例。一个小球在液体中从静止释放受到重力向下拉浮力向上推液体对颗粒的阻力也向上最终达到一个匀速沉降的终端速度。这个案例要用欧拉-拉格朗日的思路液体场作为背景连续相做简化假设我们只追踪颗粒这一条轨迹。我做了三个假设颗粒是球形且足够小颗粒间不碰撞液体是静止的。这样颗粒运动就只由重力和曳力决定。方程是牛顿第二定律mp * dv/dt (mp - mf) * g - Fdrag其中mp是颗粒质量mf是等体积液体质量等号右边第一项是净浮力重力减浮力第二项是曳力。曳力的标准形式是Fdrag 0.5 * ρf * Cd * A * v * |v|A是颗粒投影面积v是颗粒速度。这里的v*|v|不是笔误是为了让阻力方向始终与运动方向相反保证代码在颗粒反向运动时不会出现符号错误。有一个细节必须讲清楚要不要加附加质量。颗粒在液体中加速时会带动周围一部分液体一起运动这部分惯性相当于给颗粒增加了一个等效质量通常取颗粒排开液体质量的一半。稳态终端速度不受附加质量影响但它会让颗粒更慢半拍到达终端速度。模拟沉降过程时我建议加上它结果更贴近真实。3.2 参数计算从粒径到雷诺数我们算一个具体例子。石英砂颗粒直径0.5 mm密度2650 kg/m³在水里沉降。水的密度1000 kg/m³黏度0.001 Pa·s。先用斯托克斯公式估算终端速度。斯托克斯沉降公式是vt_stokes (ρp - ρf) * g * d² / (18 * μ)代入数值d是颗粒直径算出来约0.225 m/s。但这时要马上检查雷诺数Re ρf * vt * d / μ 1000 * 0.225 * 0.0005 / 0.001 ≈ 112雷诺数112已经远大于斯托克斯公式的适用范围Re远小于1所以Stokes结果不能用。这就是我强调无量纲数的原因——如果看一眼Re就不会被一个看似合理的解析值骗过去。改用Schiller-Naumann曳力修正Cd (24/Re) * (1 0.15 * Re^0.687)这个公式在Re从0.1到1000范围内都算靠谱。因为Cd和Re互相牵制需要用迭代法求终端速度。我习惯写一个很小的迭代循环先猜vt0.1 m/s迭代几次后收敛到约0.08 m/s。这个值比斯托克斯的0.225小很多这就是阻力修正的威力。3.3 完整代码与逐段解读下面这段代码是我实际在用的版本加了附加质量也做了理论终端速度的迭代计算。复制到MATLAB里直接能跑。clear; clc; close all; %% 参数区 dp 0.5e-3; % 颗粒直径m rhoP 2650; % 颗粒密度kg/m3 rhoF 1000; % 液体密度kg/m3 muF 1e-3; % 液体动力黏度Pa.s g 9.81; % 重力加速度m/s2 Vp pi * dp^3 / 6; % 颗粒体积 Ap pi * dp^2 / 4; % 颗粒投影面积 mp rhoP * Vp; % 颗粒质量 mAdd 0.5 * rhoF * Vp; % 附加质量 M mp mAdd; % 参与加速的等效质量 dt 1e-4; % 时间步长s T 0.8; % 总模拟时长s Nstep round(T / dt); %% 主循环欧拉法推进 v zeros(Nstep, 1); % 速度向下为正 y zeros(Nstep, 1); % 沉降距离 v(1) 0; y(1) 0; for k 1:Nstep-1 Re rhoF * abs(v(k)) * dp / muF; % 实时雷诺数 Cd 24 / Re * (1 0.15 * Re^0.687); % Schiller-Naumann Fdrag 0.5 * rhoF * Cd * Ap * v(k) * abs(v(k)); % 曳力符号自动与速度相反 Fbuoy (rhoP - rhoF) * Vp * g; % 净浮力 a (Fbuoy - Fdrag) / M; % 加速度 v(k1) v(k) a * dt; y(k1) y(k) v(k1) * dt; end %% 终端速度理论值迭代计算 vt 0.1; % 初猜 Cd 0; for it 1:100 Re rhoF * vt * dp / muF; Cd 24 / Re * (1 0.15 * Re^0.687); vtNew sqrt(4 * (rhoP - rhoF) * g * dp / (3 * rhoF * Cd)); if abs(vtNew - vt) 1e-8, break; end vt vtNew; end fprintf(数值终端速度%.4f m/s\n, v(end)); fprintf(理论终端速度%.4f m/s\n, vt); %% 绘图 t (0:Nstep-1) * dt; figure; subplot(2,1,1); plot(t, v, b-, LineWidth, 1.5); hold on; yline(vt, r--, LineWidth, 1.2); grid on; xlabel(时间 (s)); ylabel(速度 (m/s)); legend(模拟速度,理论终端速度,Location,southeast); title(颗粒沉降速度随时间变化); subplot(2,1,2); plot(t, y * 1000, k-, LineWidth, 1.5); grid on; xlabel(时间 (s)); ylabel(沉降距离 (mm)); title(颗粒沉降距离随时间变化);代码里最需要注意的是曳力符号处理。我写的是Fdrag 0.5 * rhoF * Cd * Ap * v(k) * abs(v(k))因为速度初始为正向下曳力向上所以加速度公式里是Fbuoy - Fdrag。如果颗粒反弹变成负速度这个写法也会自动让曳力反向正好符合物理意义。时间步长为什么取1e-4秒这个案例主要受加速度和速度变化限制。颗粒终端速度约0.08 m/s从零加速到终端速度大约需要0.1到0.2秒dt0.0001秒意味着一个时间常数内有两千个步长精度足够。再大一点比如dt0.01秒曲线会明显出现锯齿速度甚至可能发散。3.4 结果分析看到曲线后要心里有数跑完这段代码你会看到速度曲线从零迅速爬升然后慢慢平坦最终逼近理论终端速度那条红虚线。我实测的结果模拟到0.8秒时数值速度与理论值的相对偏差在1%以内说明时间步长和阻力模型都没问题。沉降距离曲线则是一个抛物线上升的形态后期接近一条直线斜率就是终端速度。这里有一个检查技巧取曲线后段两点算斜率应该和理论终端速度一致。如果不一致优先检查单位制。这个案例有很大扩展空间。把单颗粒写成多颗粒速度向量变成矩阵就能模拟几百个粒径随机分布的颗粒同时沉降。每个颗粒的速度更新互不影响MATLAB矩阵一次性就能更新。计算社会颗粒的沉降速度分布统计出来画直方图效果非常直观。很多DPM模型的第一步就是这么来的。4. 案例二垂直管内气液两相流的漂移通量模型4.1 物理背景和模型假设第二个案例换个思路回到欧拉-欧拉的工程框架解决一个实际得多的场景垂直管道内气液两相流动。垂直管里同时通入水和空气气液两相会呈现不同的流型气量小时是气泡均匀分布的气泡流气量大时小气泡聚并形成段塞流再大就变成环状流。工程上关心的是给定气液流量管道里的空隙率气相体积分数是多少总压降是多少因为这直接决定了泵怎么选、管线设计余量该留多少。漂移通量模型是处理这类问题最经典的简化模型由Zuber和Findlay在1965年提出。它的核心思想是气相实际速度不是混合物速度j那么简单而是等于分布参数C0乘以混合速度j再加上一个局部滑移速度ud。公式写出来是ug C0 * j ud这里j Jg Jl是表观气速和表观液速之和也就是混合物体积通量。C0反映的是气相在管道截面上分布不均匀带来的效应气速高时气泡会在管中心聚集跑得更快泡状流通常取1.2。ud则反映的是单个气泡受浮力作用在静止液体中的上升速度这就是漂移二字的由来。泡状流的漂移速度用Zuber-Findlay公式ud 1.53 * (σ * g * (ρl - ρg) / ρl²)^0.25σ是表面张力。这一项表面张力怎么进来了因为气泡的上升速度受气泡形状和界面效应影响表面张力起着约束气泡变形的作用。这个公式量纲上是合理的算出来后单位是m/s。有了ug之后空隙率就能直接算α Jg / ug Jg / (C0 * j ud)这是模型最漂亮的地方不需要求解微分方程一个代数公式就把空隙率确定下来了。4.2 MATLAB实现与参数选择我给的工况是管径50 mm液体表观速度约0.05 m/s气体表观速度从0.01扫到1.0 m/s研究不同气量下的空隙率和压降趋势。物性取常温常压水和空气。计算压降时分两部分重力压降和摩擦压降。重力压降很简单就是混合物密度乘重力。摩擦压降用达西公式但雷诺数和摩擦系数要用两相混合参数去算。雷诺数Re_m ρm * j * D / μm其中混合物密度ρm α*ρg (1-α)ρl混合物黏度μm αμg (1-α)*μl。摩擦系数按雷诺数分支层流取64/Re湍流用Blasius公式0.3164/Re^0.25。clear; clc; close all; %% 参数区 D 0.05; % 管径m A pi * D^2 / 4; % 截面积m2 rhoL 1000; % 液相密度kg/m3 rhoG 1.2; % 气相密度kg/m3 muL 1e-3; % 液相黏度Pa.s muG 1.8e-5; % 气相黏度Pa.s sigma 0.07; % 表面张力N/m g 9.81; % 重力加速度 C0 1.2; % 泡状流分布参数 L 10; % 管长m Ql 1e-4; % 液相体积流量m3/s Jl Ql / A; % 表观液速约0.051 m/s JgList linspace(0.01, 1.0, 60); % 表观气速扫描 % 漂移速度 ud 1.53 * (sigma * g * (rhoL - rhoG) / rhoL^2)^0.25; fprintf(漂移速度 ud %.4f m/s\n, ud); %% 主循环扫描气速 alpha zeros(size(JgList)); % 空隙率 dPdz zeros(size(JgList)); % 单位管长总压降 for k 1:numel(JgList) Jg JgList(k); j Jg Jl; % 空隙率 a Jg / (C0 * j ud); if a 0.95, a 0.95; end % 物理限制 if a 0.01, a 0.01; end alpha(k) a; % 混合物物性 rhoM a * rhoG (1 - a) * rhoL; muM a * muG (1 - a) * muL; % 雷诺数与摩擦系数 ReM rhoM * j * D / muM; if ReM 2300 f 64 / ReM; else f 0.3164 * ReM^(-0.25); end % 压降重力 摩擦 dPdz_gravity rhoM * g; dPdz_fric f * rhoM * j^2 / (2 * D); dPdz(k) dPdz_gravity dPdz_fric; end %% 绘图 figure; subplot(2,1,1); plot(JgList, alpha, b-, LineWidth, 1.5); grid on; xlabel(表观气速 Jg (m/s)); ylabel(空隙率 α); title(空隙率随表观气速变化); subplot(2,1,2); plot(JgList, dPdz, r-, LineWidth, 1.5); grid on; xlabel(表观气速 Jg (m/s)); ylabel(单位管长压降 (Pa/m)); title(压降随表观气速变化);运行这段代码之前先手动算一下漂移速度的量级。把σ0.07、g9.81、密度差代入算出来ud大约是0.2 m/s左右和文献里气泡在水中的终端速度是一个量级这就合理。如果算出来ud等于2 m/s或者0.02 m/s那一定是单位制出了问题。4.3 结果解读为什么压降曲线有个坑跑完代码你会看到两个图。第一个图α随Jg增大而增大但斜率逐渐变小。当Jg很大时α趋向于1/1.2约0.83因为模型里C01.2气速再大空隙率也被锁在这个上限附近。这对应实际物理情形段塞流里液相不可能被完全排干总有液膜和液塞存在。第二个图更有意思。压降并不是单调的而是先下降降到一个极小值后再上升。这是垂直气液两相管流的经典现象小气量刚开始增加时原本静止的液相被气流撬动空隙率增加混合物密度降低重力压降减少的效应超过了摩擦增加所以总压降下降。到某个气速之后气相速度过高壁面摩擦压降开始占主导总压降掉头向上。工程上这个压降极小值对应的点是多相流动不稳定性的根源。气举工艺设计时如果操作点落在这个区域很容易出现液塞-气速耦合振荡产量忽高忽低。看到这个坑你就知道这个几十行的MATLAB模型已经在帮你抓住工程问题的关键趋势了。4.4 这个模型的边界在哪里必须诚实地说明这个模型只适用于泡状流和较低空隙率区间。我代码里设了α上限0.95就是为了避免在极高气速下模型失真。气速超过一定值后流型转入环状流C0和ud都不再是常数甚至要换一套经验关联式。MATLAB的优势在这里就体现出来了想修模型直接在参数区改个公式再跑一遍十秒钟看趋势变化。另外摩擦压降用的是简化单相等效方法严格意义上应该用Lockhart-Martinelli两相摩擦乘子来处理。但对于教学和趋势分析这个简化已经能反映出最核心的物理特征。真要做工程设计用这个程序定初始参数再用专业软件复核效率会高很多。5. 常见问题与排查技巧实录代码跑不起来或者结果离谱通常不是什么玄学问题就是下面几类。5.1 颗粒速度发散发散的第一嫌疑是时间步长过大。颗粒模拟里能量方程是显式推进dt太大就数值不稳定。判断依据很简单把dt减半如果结果变化很大说明时间步长没收敛。第二个嫌疑是曳力计算时速度出现零值导致Re0然后Cd24/0直接爆掉。我在代码里用abs(v(k))但如果速度恰好从零开始第一步Re0就会出问题。稳妥的做法是给Re加一个小值保护Re rhoF * abs(v(k)) * dp / muF 1e-30。5.2 模拟结果与理论值对不上先别急着怀疑代码先自己手算一遍。颗粒沉降案例里理论终端速度用迭代法能算出来你把数值结果和理论值分别打印出来对比如果偏差超过5%百分之九十九是阻力模型或者单位制的问题。我见过最典型的错误是把直径0.5 mm直接当成0.5代入公式算出来终端速度大了1000倍。5.3 压降曲线出现负值或跳变压降突然跳变通常发生在流型转换区域。如果摩擦系数在层流和湍流分支之间来回切换曲线会出现锯齿。这是物理模型本身的问题不是程序bug可以加一个平滑过渡或者在参数区限制Re范围。另外空隙率α设置了0.95的上限后压降在接近上限时会明显失真论文里画图时要标注适用范围。5.4 MATLAB计算很慢多相流模拟的慢分为两种。一种是真慢比如用双层循环处理大网格每步都重新分配变量这种慢可以通过预分配数组、to矩阵操作来根治。另一种是本来就该慢你非要拿MATLAB跑一个三维VOF网格那不属于优化能解决的范畴。遇到后者及时换工具。5.5 画动画卡成幻灯片用drawnow逐帧更新会导致动画非常卡。我用的是drawnow limitrate它能自动丢帧保证界面响应。还有一个技巧先把计算结果全部存储到预分配的三维数组里仿真结束后统一绘图而不是边算边画这样既能保证计算速度又能随时转出任意时刻的云图。我把常见问题整理成一个速查表方便以后排查现象可能原因排查顺序速度曲线发散时间步长过大Re为0减小dt检查初始Re终端速度偏大阻力模型用错或单位制错手算理论值检查单位空隙率超过1气相速度公式或表观气速错误检查C0和ud量级压降随气速无变化摩擦系数没切换层湍流分支打印ReM检查分支逻辑计算慢循环没有向量化预分配改矩阵运算动画卡drawnow刷帧太频繁用drawnow limitrate6. 让我最后再啰嗦几句我自己的体会是写多相流程序时最花时间的地方从来不是写代码而是把每个物理量为什么要这样定义、每个参数为什么要取这个值想清楚。代码只是把物理模型翻译成数值操作而已。MATLAB的好处是翻译成本很低改一个阻力系数公式重跑一遍就是几秒钟的事这种快速试错的过程能帮助你形成很多直觉。最后分享一个我习惯用的工作流拿到一个多相流问题先不出手写代码在纸上把连续相、分散相、主导力、无量纲数列出来然后用最简单的手算公式估算一个预期量级。比如颗粒沉降先算终端速度管道流先算空隙率大致范围。手算和程序结果对上再往下扩展模型。这样你永远知道自己算出来的结果倒不靠谱也算是对自己时间负责。