做流体数值计算的人应该都经历过这种尴尬速度云图拉出来整片颜色接近哪里是涡、哪里是剪切层、涡核中心在哪、旋转强弱如何光靠盯屏幕根本判断不了。真正能让我们一眼读懂流场结构的反而是涡量场。这篇博文用MATLAB把通道流中瞬时涡量场的算法实现完整走一遍从数学定义、数值差分到可视化出图逐段拆源码项目代码能直接拿去复现也能作为PIV后处理、CFD结果诊断的通用模板。适合正在学流体数值计算、做MATLAB可视化或者想快速验证自己流场数据的人读完可以直接抄作业。1. 瞬时涡量场先搞清楚我们算的是什么1.1 涡量的物理含义与数学定义涡量vorticity的严格定义是速度场的旋度记作 \omega \nabla \times \mathbf{V}。在二维流动中只有垂直于流动平面的分量存在表达式可以简化为[ \omega_z \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} ]其中 u 和 v 分别是 x、y 方向的速度分量。注意这个公式的符号约定至关重要——正涡量对应数学正方向上的逆时针旋转负涡量对应顺时针旋转。很多教材里因为坐标系画法不同正负号往往会反过来初学者最容易在这里翻车。生活化地理解涡量速度场告诉你“每一片草叶往哪个方向跑”涡量告诉你“这片草叶自己转不转、转得多快”。两个流场可能速度分布完全不同但只要绕某一点的旋转强度相同涡量就接近。这也是为什么在剪切层、边界层、涡脱落这些场景里涡量云图比速度云图更具辨识度。瞬时涡量场instantaneous vorticity field强调的则是“某一时刻”的空间分布。流动研究中时均涡量反映的是平均剪切和平均旋转而瞬时涡量保留了涡结构演化过程中的细节比如涡核的拉伸、破裂、合并这些现象在时均场里会被平滑掉。本文要实现的就是针对通道流模型取某一个确定时刻的流场数据计算并可视化该时刻的涡量分布。1.2 通道流模型与坐标约定通道流channel flow也叫平面泊肃叶流是流体力学最经典的基准算例之一两平行平板之间充满流体沿 x 方向施加压力梯度驱动流动z 方向无限展向所以可以简化为二维问题。流场的主要特征是无量纲速度剖面近似抛物线分布壁面处速度为零中心线速度最大同时壁面附近存在强剪切层而剪切层正是高涡量集中的区域。坐标约定上我把 x 设为流向水平方向y 设为法向垂直于壁面方向通道上下壁面分别位于 y -1 和 y 1通道高度 H 2流向计算域长度取 Lx 4。所有量采用无量纲形式中心速度设为 Umax 1。这样做的好处是计算结果干净、单位统一也方便直接和文献中的无量纲结果对照。之所以选通道流而不是更简单的自由剪切流是因为通道流同时包含了壁面边界、中心主流、剪切层分布这几个关键要素能够把“壁面涡量生成”和“中心区涡量输运”这一对物理过程同时展示出来对理解涡量场非常有帮助。而瞬时性的体现则通过在背景剖线上叠加若干个位置随时间演化的Oseen涡来实现这样每个时刻都会得到一幅结构清晰、有变化的瞬时涡量场。2. 算法设计思路从流场到涡量场2.1 整体处理管线瞬时涡量场计算并不复杂核心思路可以拆成几条明确的处理链先构造计算网格然后生成或读取该时刻的速度场再用数值微分计算速度梯度最终得到涡量场并做可视化。完整管线如下表环节输入输出关键工具网格构造计算域尺寸、网格数网格坐标矩阵 X、Ymeshgrid流场生成涡参数、背景剖线u(x,y,t)、v(x,y,t)解析函数叠加数值微分速度场矩阵、步长\partial v/\partial x、\partial u/\partial y差分模板涡量合成两个梯度场涡量场 \omega(x,y,t)梯度相减可视化涡量场、速度场云图、矢量图、动画contourf、quiver、VideoWriter这里我把“流场生成”和“涡量计算”分开而不是直接写一个函数生成涡量场是为了模块化以后如果用户手里有真实CFD结果或PIV实验数据只需要替换第二步的数据来源后面的差分、可视化代码可以原封不动复用。这也是代码设计的核心原则之一。2.2 数值微分方案选型计算涡量的核心是对速度场求空间偏导数。MATLAB里有内置的梯度函数 gradient默认情况下它也是用差分实现的但边界点采用单边差分会带来一定误差。如果想完全掌控算法细节建议手写差分模板。以 \partial v/\partial x 为例二阶中心差分公式为[ \frac{\partial v}{\partial x}\bigg|{i,j} \approx \frac{v{i,j1} - v_{i,j-1}}{2\Delta x} ]中心差分的截断误差是 O(\Delta x^2)相比向前差分 O(\Delta x) 精度高一阶而且对正弦波的相位误差更小。对于光滑的解析流场例如本文用的叠加Oseen涡中心差分几乎能达到机器精度级别的近似对于实验测量数据中心差分还能起到轻微平滑作用不至于像向前差分那样局域偏差过大。边界处无法使用中心差分只能退而求其次采用一阶单边差分[ \frac{\partial v}{\partial x}\bigg|{i,1} \approx \frac{v{i,2} - v_{i,1}}{\Delta x} ]这一处理虽然会降低边界精度但对于通道流而言壁面本身是强剪切区域涡量值本身就很大局部梯度误差并不会严重影响整体云图结构。我在代码中单独写了边界分支让逻辑一目了然。2.3 合成流场数据为什么不用现成CFD结果很多读者会问为什么不用FLUENT或者OpenFOAM算一个通道流再导入MATLAB我的回答是对于这篇博文的目标场景合成解析流场有两个不可替代的优势。第一可复现性。解析流场不依赖任何商业软件版本、网格划分方式或迭代收敛水平只要参数固定任何机器上跑出来的结果都一致非常适合用来讲算法和调试代码。第二可控性。我想让涡心在哪个位置、强度多大、正负旋转方向如何都可以精确设定方便验证可视化结果是否符合理论预期——比如让两个正负涡量交替排列马上就能观察云图的正负色分布。合成流场的方案是“背景抛物线剖线 叠加Oseen涡”。Oseen涡也叫Lamb-Oseen涡是层流涡的解析解切向速度分布为[ u_\theta(r) \frac{\Gamma}{2\pi r}\left(1 - e^{-r^2 / (2\sigma^2)}\right) ]其中 \Gamma 是涡强度环路积分量\sigma 是涡核半径。这个速度分布在 r0 附近为零向外先增大后衰减形状非常接近真实流动中的集中涡结构。把极坐标速度转换到笛卡尔坐标时在涡心 (x_c, y_c) 处需要避免 r0 的奇点代码中我采用了 max(r2, 1e-10) 的截断处理。3. MATLAB源码实现逐段拆解3.1 主脚本骨架与参数设置先给出完整的主脚本结构。建议把文件命名为 instantaneous_vorticity_channel.m所有参数集中在文件头部方便后续修改。clear; clc; close all; %% 参数设置 Lx 4.0; % 流向计算域长度 Ly 2.0; % 法向通道高度 nx 201; % 流向网格数 ny 101; % 法向网格数 dx Lx / (nx - 1); % 流向网格间距 dy Ly / (ny - 1); % 法向网格间距 % 网格坐标X 是 ny*nxY 是 ny*nx x linspace(0, Lx, nx); y linspace(-Ly/2, Ly/2, ny); [X, Y] meshgrid(x, y); % 背景抛物线剖线参数 Umax 1.0; % 通道中心最大速度 h Ly / 2; % 半通道高度 % Oseen涡参数可自行增减数量 gamma [0.6, -0.4, 0.3]; % 涡强度正为逆时针、负为顺时针 xc0 [0.8, 1.8, 2.9]; % 涡心初始流向位置 yc0 [-0.2, 0.25, -0.15]; % 涡心法向位置 sigma [0.15, 0.12, 0.1]; % 涡核尺度 nc length(gamma); % 涡数量这里的网格分辨率选择需要解释一下。201×101 的网格对应流向网格间距约0.02法向约0.02比最小涡核尺度0.1 小约5倍。也就是说一个涡核直径范围内大概分布了10个网格点足以解析涡量梯度。如果网格太疏比如只有51×26涡核会呈现明显的锯齿状差分精度也会下降。3.2 速度场生成背景剖线叠加Oseen涡速度场的生成逻辑分为两步先给抛物线背景剖线再逐涡叠加扰动速度。%% 速度场生成 % 背景通道流抛物线速度剖线 u Umax * (1 - (Y / h).^2); v zeros(size(u)); % 叠加 Oseen 涡 for k 1:nc rx X - xc0(k); ry Y - yc0(k); r2 rx.^2 ry.^2; r2 max(r2, 1e-10); % 防止奇点 factor gamma(k) / (2*pi) .* (1 - exp(-r2 / (2*sigma(k)^2))); % Oseen涡在笛卡尔系下的速度分量 u u - factor .* ry ./ r2; v v factor .* rx ./ r2; end逐行解释这段代码。Umax * (1 - (Y/h).^2) 直接实现了无量纲抛物线剖面壁面 Y ±1 处速度为零中心 Y 0 处速度为 Umax。随后进入涡叠加循环rx、ry 是当前网格点相对涡心的距离向量。r2 max(r2, 1e-10) 的作用是避免涡心处除以零同时把 1e-10 作为极小半径下的截断。代码里的 factor 变量是整个叠加的核心。当 r2 很小时1 - exp(-r2/(2σ^2)) 近似为 r2/(2σ^2)factor ≈ Γ*r2/(4πσ^2)所以 u、v 的修正量都趋近于零涡心处速度连续当 r2 很大时指数项趋近于零factor ≈ Γ/(2π)速度衰减趋势符合远场涡的 1/r 特征。整段代码不需要循环遍历网格点而是用矩阵运算一次完成所有网格的速度修正这也是MATLAB推荐的向量化写法。3.3 涡量计算与边界差分处理速度场生成之后进入核心的涡量计算环节。这里我采用手动差分模板保证边界处理完全可控。%% 涡量计算 \omega dv/dx - du/dy omega zeros(ny, nx); % 内部区域二阶中心差分 for i 2 : ny - 1 for j 2 : nx - 1 dvdx (v(i, j1) - v(i, j-1)) / (2*dx); dudy (u(i1, j) - u(i-1, j)) / (2*dy); omega(i, j) dvdx - dudy; end end % 边界一阶单边差分 % 左右边界x 方向 for i 2 : ny - 1 omega(i, 1) (v(i, 2) - v(i, 1)) / dx - (u(i1, 1) - u(i-1, 1)) / (2*dy); omega(i, nx) (v(i, nx) - v(i, nx-1)) / dx - (u(i1, nx) - u(i-1, nx)) / (2*dy); end % 上下边界y 方向 for j 2 : nx - 1 omega(1, j) (v(2, j) - v(1, j)) / (2*dx) - (u(2, j) - u(1, j)) / dy; omega(ny, j) (v(ny, j) - v(ny-1, j)) / (2*dx) - (u(ny, j) - u(ny-1, j)) / dy; end % 四个角点直接用邻近值填充避免出现除零 omega(1, 1) omega(2, 2); omega(1, nx) omega(2, nx-1); omega(ny, 1) omega(ny-1, 2); omega(ny, nx) omega(ny-1, nx-1);这里我需要强调两个重点。第一MATLAB 中 meshgrid(X, y) 生成的矩阵第一维是 y 方向行数 ny第二维是 x 方向列数 nx所以 u(i1, j) 对应法向相邻点u(i, j1) 对应流向相邻点。这一点如果搞反算出来的涡量会整体错位。第二差分顺序。\partial v/\partial x 中 v 是 y 方向速度分量需要对 x 求偏导因此索引变化在第二维\partial u/\partial y 中 u 是 x 方向速度分量需要对 y 求偏导因此索引变化在第一维。代码中我严格按照这个约定实现这也是最容易写错的地方。边界处理上由于中心差分需要两侧点在左边界 j1 处改用前向差分右边界 jnx 处改用后向差分上下边界同理。角点缺少两个方向的差分信息直接用邻近内部点的值填充。对于瞬时涡量可视化而言角点区域不是关注重点这一简化完全足够。3.4 让瞬时场真正“瞬时”加入时间演化上面的代码已经能算出一个确定时刻的瞬时涡量场但既然题目强调“瞬时”我建议把时间维加进来观察涡量场随时间的演化。思路不复杂让涡心位置以某个速度做对流运动同时涡强度可以做小幅震荡每个时间步重新生成速度场并计算涡量。%% 瞬时演化多个时刻循环 tlist linspace(0, 4, 81); % 时间序列 Uc 0.6; % 涡心整体对流速度 A_osc 0.15; % 涡强度振荡幅值 figure(Color, w, Position, [100 100 900 420]); for it 1 : length(tlist) t tlist(it); % 更新涡心位置和强度 xc xc0 Uc * t; gam gamma .* (1 A_osc * sin(2*pi*t (1:nc))); % 重新生成速度场 u Umax * (1 - (Y / h).^2); v zeros(size(u)); for k 1 : nc rx X - xc(k); ry Y - yc0(k); r2 max(rx.^2 ry.^2, 1e-10); factor gam(k) / (2*pi) .* (1 - exp(-r2 / (2*sigma(k)^2))); u u - factor .* ry ./ r2; v v factor .* rx ./ r2; end % 计算涡量复用前面的差分模块 omega zeros(ny, nx); for i 2 : ny-1 for j 2 : nx-1 omega(i, j) (v(i, j1) - v(i, j-1)) / (2*dx) ... - (u(i1, j) - u(i-1, j)) / (2*dy); end end % 绘图与刷新 clf; contourf(X, Y, omega, 40, LineColor, none); hold on; hq quiver(X(1:5:end, 1:5:end), Y(1:5:end, 1:5:end), ... u(1:5:end, 1:5:end), v(1:5:end, 1:5:end), 1.2); hq.Color [0.2 0.2 0.2]; colormap(blue_white_red(256)); caxis([-1 1]); axis equal; xlim([0 Lx]); ylim([-h h]); title(sprintf(Instantaneous Vorticity Field, t %.2f, t)); xlabel(x); ylabel(y); colorbar; drawnow; % 此处可保存视频帧见第4节 end时间演化这里有个物理合理性需要说明。真实通道流中涡的对流速度通常不是均匀的——中心区快、近壁区慢甚至存在逆向运动。我这里统一用 Uc 0.6 做整体平流是简化处理目的是演示算法框架。如果接入真实流场数据时间演化逻辑完全不需要改只需要把“重新生成速度场”替换为“读取下一时刻数据”。4. 可视化关键环节怎样把瞬时涡量画得科学又直观4.1 绘图函数选型contourf、pcolor 还是 surfMATLAB 里画标量场云图常见的有 pcolor、contourf、surf、imagesc 四种。很多人直接拿 pcolor 画结果发现图片边缘有网格线、颜色过渡不连续其实是用错了场景。函数适用场景优点缺点contourf标量场填充云图平滑、可控制等值线数量默认去网格线需要指定填充等级数pcolor标量场伪彩图原生网格化显示边缘有网格线需 shading interp 处理surf三维曲面加高度显示适合展示梯度起伏二维俯视时优势不明显imagesc均匀网格标量图速度快、内存小坐标缩放需额外处理我的建议是静态云图优先用 contourf配合 LineColor 设为 none可以得到干净的填充图。如果想让涡量场显示立体感再用 surf 加 view(2) 视角。本文演示采用 contourf quiver 的组合云图展示涡量大小和方向矢量箭头展示当地速度方向二者叠加能让读者直观感受“涡量集中在剪切层”这一物理事实。4.2 颜色映射不要默认 jet用蓝白红发散色涡量有正有负如果用默认 jet 色图零涡量区域会被映射到绿色正负涡量无法形成直观的视觉对比。正确的做法是采用蓝-白-红发散色图负涡量蓝色、零涡量白色、正涡量红色。这样流场中顺时针涡、逆时针涡、无旋区域一眼就能分辨。function cmap blue_white_red(n) % 生成蓝-白-红发散colormap零值居中 if nargin 1, n 256; end top [0 0 0.6; 0 0 1; 1 1 1; 1 0 0; 0.6 0 0]; cmap interp1(linspace(0, 1, size(top, 1)), top, linspace(0, 1, n)); end这段自定义函数用 interp1 对五个控制点做线性插值生成 256 色的渐变映射。控制点采用深蓝-纯蓝-白色-纯红-深红的结构这样色图两端是深色、中间是浅色动态范围均匀。需要特别指出的是 caxis 的取值如果直接用默认范围某个极强涡量值会把整个色带拉伸导致弱涡结构完全看不清。建议设定对称范围比如 caxis([-1 1])再根据实际涡量幅值调整。4.3 图面精细控制与出图出图阶段有几个小细节经常被忽略。第一colorbar 一定要加标签标明“vorticity ω”并注明无量纲单位。第二quiver 的箭头要适当抽稀否则密集网格下矢量箭头会糊成一片。我在代码里用 1:5:end 抽稀每隔 5 个网格点画一个箭头速度缩放因子设为 1.2这个参数可以根据视野大小微调。第三坐标轴纵横比要设置 axis equal避免通道流看起来被拉扁或拉伸。最后用 exportgraphics 输出高分辨率图片exportgraphics(gcf, instantaneous_vorticity.png, Resolution, 300);提示exportgraphics 在 MATLAB R2020a 及以上版本可用老版本可以使用 print(gcf, -dpng, -r300, instantaneous_vorticity.png)效果几乎相同。5. 常见问题与排查技巧实录5.1 典型问题速查表以下问题是我在实际编码和帮人调试过程中遇到频率最高的整理成速查表现象可能原因排查方法涡量云图全为 0速度场 u、v 写反涡量计算索引错误打印 dvdx 和 dudy 检查是否非零涡量颜色只有单一色colormap 选错或 caxis 范围过大修改为蓝白红发散色caxis 设为涡量幅值对称范围涡核出现明显锯齿网格分辨率不足增加 nx、ny使涡核尺度覆盖 5 个以上网格点边界涡量异常高差分精度不足或流场本身物理合理用解析涡验证边界是否弱信号而非计算错误图像边缘大量黑线contourf 未设 LineColor none加上 LineColor 设置为 nonequiver 箭头方向与涡旋方向不符速度差分符号约定错误或者涡强度 \Gamma 正负定义与预期不一致用单个Oseen涡验证可视化方向动画过程中涡心跳出计算域时间推进距离过大调整 tlist 范围或增大计算域 Lx颜色图对比度低等值线数量太少或数据本身动态范围小增加 contourf 的等值线数量到 40 以上5.2 三个我踩过的坑第一个坑是速度分量和索引维度不匹配。早期我写代码时用 size(u) 确认维度后以为 u(i, j1) 是 x 方向增量结果在 meshgrid 约定下第二维才是 x第一维是 y导致涡量场整体错位。排查了很久最后是把单个 Oseen 涡的解析解和数值计算比对才定位到问题。建议以后拿到网格后先打印 X(1:3, 1:3) 和 Y(1:3, 1:3) 观察维度方向再开始写差分。第二个坑是 caxis 的动态范围。当流场中存在较强剪切层时壁面附近的涡量值往往比中心区域的涡高出几十倍如果直接用默认 colorbar 范围中心区域的弱涡会被“冲淡”到一片白色看不出结构。解决方法是把 caxis 设为对称区间或者用百分位截断比如设定 caxis([-prctile(abs(omega(:)), 95), prctile(abs(omega(:)), 95)])只显示95%范围内的数据避免离群值主导色带。第三个坑是正负涡量的符号判断。流体力学中逆时针旋转是否对应正涡量完全取决于坐标系的定义和速度分量的符号约定。如果不小心把 v 对 x 的偏导减反了整个云图的正负区域会完全对调表面上看起来像“镜像”但物理意义完全不同。验证是否算反的方法很简单单独放一个强度为正的Oseen涡如果涡量云图的中心确实显示正红色反向则说明符号反了。6. 瞬时涡量场的进一步应用6.1 时间序列动画与定量统计本文第3.4节的循环已经展示了瞬时涡量场的时间演化但严格来说动画还只是定性观察。想要定量分析涡量场的演化规律可以从时间序列中提取几个关键统计量。第一个是全域总涡量强度定义为 \int |\omega| dA反映整个区域旋转强度的总和。计算时可以用矩阵求和乘以网格面积近似sum(abs(omega(:))) * dx * dy。如果这个值随时间持续增长说明流动中存在持续的旋涡增强机制如果衰减说明黏性耗散占主导。第二个是涡量极值位置追踪。在每个时间步找到 |\omega| 最大值对应的坐标记录其轨迹。这个轨迹可以用来间接反映涡心的运动路径尤其适用于识别涡的合并或分裂事件。实现上利用 find 函数定位即可不需要复杂算法。6.2 从CFD或PIV数据接入的改法真实工程场景里速度场往往来自Fluent、OpenFOAM、STAR-CCM等CFD软件或者PIV粒子图像测速实验。这些数据导入MATLAB后核心涡量计算函数完全不需要改只需要处理数据接口。以从CFD导出的CSV文件为例表格包含列x, y, u, v接入逻辑如下% 假设已有 data.mat包含 x、y、u、v 四个列向量 data load(cfd_instant.csv); xcol data(:, 1); ycol data(:, 2); ucol data(:, 3); vcol data(:, 4); % 将散点数据 reshape 为网格矩阵 nx length(unique(xcol)); ny length(unique(ycol)); Xg reshape(xcol, ny, nx); Yg reshape(ycol, ny, nx); U reshape(ucol, ny, nx); V reshape(vcol, ny, nx);随后直接把 U、V 传给涡量差分模块即可。这里有个前提CFD导出的数据必须按 y 排列优先先固定一行 y再循环 x否则 reshape 会错位。检查方法是在 reshape 前用 unique 看 x、y 的排列顺序。注意如果CFD网格是非均匀的二阶中心差分公式需要改为非均匀网格形式\partial u/\partial y 的差分系数不再是 1/(2dy)而是由局部网格间距决定。这个情况会复杂一些建议用MATLAB的 gradient 函数先跑通流程再考虑高精度替代方案。6.3 从实验数据到涡量场的几个额外提醒实验测量数据尤其PIV与合成或CFD数据有本质区别噪声水平高、空间分辨率有限、可能存在不少错误矢量。直接用本文的中心差分处理会放大噪声导致涡量云图出现“盐粒状”噪点。我的处理经验是在计算涡量前先对速度场做一次适度的空间平滑常见做法是中值滤波medfilt2或高斯滤波imgaussfilt滤波核大小建议不超过5×5否则会过度抹平小尺度涡结构。平滑之后再套用本文的差分代码得到的涡量场会干净得多。代价是小尺度涡量峰值的幅值会被压低所以如果做定量分析需要在文中明确说明平滑对结果的影响。我自己处理PIV数据的固定流程是原始速度场 - 离群矢量剔除基于中值检验 - 3×3高斯平滑 - 中心差分 - 涡量计算 - 发散色云图出图。这套流程处理室内小尺度水槽实验数据已经跑了大半年效果稳定。结尾做流场可视化这几年我最大的体会是涡量不是速度场的附属品而是真正能揭示流动本质的物理量。速度场给人的是“运动感”涡量场给人的是“结构感”。一个瞬时流场一旦换成涡量视角旋涡的生成、迁移、耗散过程就全都清晰了。最后分享一个小经验任何时候拿到一个新的流场数据处理任务我都建议先跑一次“合成解析场验证”——用一个已知解析解比如本文的Oseen涡测试自己的差分和可视化代码确认结果与理论完全一致后再处理真实数据。这一步往往能帮你省下后面几天调试时间。代码部分可以直接复制到MATLAB里运行参数改动也都在文件头部大家可以根据自己的流场尺寸做适配。