企业官网建设流程全解析

做流体实验或者跟CFD打交道的人，第一次去画瞬时涡量场基本都会翻车。不是算不出来，而是画出来跟论文里的涡量云图完全两个物种：要么噪点炸裂像一地芝麻，要么等值线密集到看不出任何结构，更别提上下壁面附近的边界条处理，一没留意就是一条条横贯全场的伪影。我自己在MATLAB里从PIV数据、DNS输出到解析流场都折腾过不少轮，这篇文章就把平面通道流（plane channel flow）里瞬时涡量场的算法实现和MATLAB源码一次讲透。文章核心关键词只有一个：瞬时涡量场，但会牵出数值差分、边界格式、可视化参数这些躲不开的配套问题。适合流体方向的研究生、做PIV实验的工程师，以及刚入门CFD可视化、想在MATLAB里快速复现涡结构云图的人。

1. 瞬时涡量场到底在算什么

1.1 涡量的定义没那么玄

涡量就是速度场的旋度，数学上写作 ω = ∇ × u。这是流体力学里衡量流体微团旋转强度的核心量，跟“涡旋”是两回事，不需要流线转圈才有涡量，直剪切流照样有涡量。

在二维平面通道流里，速度场写成 u = (u(x, y), v(x, y), 0) ，这时候涡量只剩一个非零分量，也就是垂直纸面的分量：

ω_z = ∂v/∂x − ∂u/∂y

这个公式是所有后续算法的起点。注意这里有两个偏导数，一个对x求，一个对y求，二者相减，方向性决定了正负涡量代表相反旋转方向。

瞬时涡量场，字面意思就是某个时刻 t 对应的涡量空间分布 ω_z (x, y, t)。因为是瞬时态，场里通常同时存在多个尺度的结构：大尺度剪切层贡献的、边界层内高频小涡贡献的、局部扰动产生的，这些都是瞬时场的一部分，和后处理做时间平均得到的时均涡量场完全不同。

1.2 为什么通道流特别适合拿来讲瞬时涡量

平面通道流是经典壁面剪切流动，两块无限大平板中间夹着流体，主流方向为x，法向为y。这种流场看似简单，却是理解转捩、湍流拟序结构、近壁条带结构的最佳简化模型。瞬时涡量场能直接暴露流场中的“涡结构”，比单纯看速度矢量场要直观得多。

很多人在这一步会有一个误解：直接用速度场画quiver，看着箭头绕圈就认为涡在场中。其实旋转箭头密集的地方当然是涡量大的地方，但通道流基流本身也有很大的涡量——抛物线速度沿法向的梯度很大，处处不为零。于是画出来的图里，中心区颜色淡、近壁区两道深色高涡量带，局部真实的小涡结构反而被掩盖掉了。这就要靠算法区分“基流贡献的涡量”和“扰动涡量”，处理方法后面源码部分会细讲。

2. 数值差分的格式选择与边界处理

2.1 中心差分是最稳妥的默认答案

如果手里有一组已知速度场数据，不管是PIV测的还是DNS导出的，涡量计算第一步都是用有限差分近似两个偏导数。

对等距网格，最常用的二阶中心差分公式是：

(∂u/∂y){i,j} ≈ (u{i,j+1} − u_{i,j−1}) / (2Δy)

(∂v/∂x){i,j} ≈ (v{i+1,j} − v_{i−1,j}) / (2Δx)

其中 (i, j) 是网格点下标，Δx、Δy是网格步长。这套格式的截断误差是二阶量级，在绝大多数CFD后处理和PIV后处理场景里足够用。它最大的优点是无偏、能保符号对称性，不会系统性放大某一方向的数值假象。

为什么默认用中心差分而不是向前或向后差分？因为前向差分的精度只有一阶，而且会引入方向的系统偏差，在一个对称的通道流里，这种偏差会导致涡量场上下不对称，这在后续分析里会非常恼人。我测试过头，相同数据用前向差分算出来的近壁峰值，位置跟真实值能偏差好几个百分点。

2.2 边界点不能用中心差分，必须单边格式

通道流有上下两个壁面，中心差分在边界点 i, j=1 处需要用到 j=0 的节点，这个节点不存在。直接在两侧各削掉一行数据，很多人都会这么干，但削掉之后近壁涡量峰值就丢了，而峰值恰恰是壁面剪切的直接体现。

正确做法是在上下边界使用二阶单边差分，比如在 y=0 处：

(∂u/∂y){i,1} ≈ (−3u{i,1} + 4u_{i,2} − u_{i,3}) / (2Δy)

这就是二阶精度单边差分，原理是构造一个过三点 (y1, y2, y3) 的二次插值多项式，再求导。在 y=ymax 处对称处理即可。这样近壁涡量算出来精度跟内部点一致，不会出现边界断层。

2.3 网格分辨率不足会直接漏掉小尺度涡

这一步最容易被忽略。涡量包含了流场的速度梯度信息，对小尺度结构天然敏感。如果网格太粗，两个相邻网格点之间夹着一个直径只有一个网格宽度的小涡，差分结果会把这个小涡平滑得几乎看不见。反之，如果测量数据本身含高频噪声，高分辨率差分反而会放大噪声，造成涡量云图满屏“椒盐”。

经验法则：要可靠分辨尺度为 L_s 的涡结构，网格尺寸至少要小于 L_s/3。PIV数据通常用互相关窗口大小作为有效分辨率，别被采集间隔骗了。

下面这个表是我处理不同类型数据时用的默认方案：

3. MATLAB源码逐段拆解

3.1 用解析流场做算法验证

我不建议直接拿真实测量数据上来就调可视化参数，因为真实数据的“真值涡量”是未知的，出问题很难判断是算法错了还是数据本身的问题。更好的路径是先用解析构造的流场走通全流程，验证算法精度，再换真实数据。

这套示例的流场构造思路是：一个抛物线基流叠加两个高斯型局部涡扰动。基流模拟通道流层流解，高斯涡模拟瞬时流场里的局部涡结构。这样瞬时涡量场的解析表达式是已知的，可以直接对比数值计算结果。

% 构造平面通道流的瞬时速度场 % 网格参数 Lx = 2*pi; % 流向长度 H = 2; % 通道高度（y从-H/2到H/2） nx = 128; ny = 128; x = linspace(0, Lx, nx); y = linspace(-H/2, H/2, ny); [X, Y] = meshgrid(x, y); % 基流：抛物线泊肃叶解 U0 = 1.0; u_base = U0 * (1 - (2*Y/H).^2); v_base = zeros(size(X)); % 扰动流函数：两个高斯涡 A = 0.15; % 扰动量级 sigma1 = 0.12; sigma2 = 0.08; psi1 = A * exp(-((X - Lx/3).^2 + (Y - H/6).^2) / (2*sigma1^2)); psi2 = -A * exp(-((X - 2*Lx/3).^2 + (Y + H/6).^2) / (2*sigma2^2)); psi_pert = psi1 + psi2; % 扰动速度由流函数求导得到 [~, psi_y] = gradient(psi_pert, x(2)-x(1), y(2)-y(1)); [psi_x, ~] = gradient(psi_pert, x(2)-x(1), y(2)-y(1)); u = u_base + psi_y; % u = du_base/dy? 不，这里psi_y = dpsi/dy v = v_base - psi_x; % v = -dpsi/dx

这里用gradient求流函数的偏导，得到扰动速度后叠加到基流，就得到了一个带局部涡结构的瞬时通道流速度场。解析涡量为基流涡量加扰动涡量，扰动涡量由 ω_pert = −∇²ψ 解析给出：

ω_pert 的理论表达式需要对高斯函数二次求导，示例里直接用数值拉普拉斯或者等差分得到近似，然后对比验证。

3.2 核心函数：二维涡量计算

下面这个函数就是全篇的核心工具，输入二维速度场和网格步长，输出瞬时涡量场。内部做了内部点中心差分、边界点单边差分，同时支持可选的前置平滑。

function omega = vorticity2d(u, v, dx, dy) % 计算二维流场的涡量场 omega_z = dv/dx - du/dy % 输入：u,v - 二维速度分量矩阵（尺寸一致），dx,dy - 网格步长 % 输出：omega - 涡量场 [ny, nx] = size(u); omega = zeros(ny, nx); % 内部点：二阶中心差分 omega(2:end-1, 2:end-1) = ... (v(2:end-1, 3:end) - v(2:end-1, 1:end-2)) / (2*dx) - ... (u(3:end, 2:end-1) - u(1:end-2, 2:end-1)) / (2*dy); % 上下壁面（y方向边界）：二阶单边差分 % du/dy 在 y=1 处：(-3u1 + 4u2 - u3)/(2dy) omega(1, 2:end-1) = ... (v(1, 3:end) - v(1, 1:end-2)) / (2*dx) - ... (-3*u(1, 2:end-1) + 4*u(2, 2:end-1) - u(3, 2:end-1)) / (2*dy); omega(end, 2:end-1) = ... (v(end, 3:end) - v(end, 1:end-2)) / (2*dx) - ... (3*u(end, 2:end-1) - 4*u(end-1, 2:end-1) + u(end-2, 2:end-1)) / (2*dy); % 左右入口出口边界（x方向）：这里假设周期边界或简单外插 % 先用零梯度外插，保证整图没有空白 omega(2:end-1, 1) = omega(2:end-1, 2); omega(2:end-1, end) = omega(2:end-1, end-1); % 四个角平均处理 omega(1,1) = (omega(2,1) + omega(1,2)) / 2; omega(1,end) = (omega(2,end) + omega(1,end-1)) / 2; omega(end,1) = (omega(end-1,1) + omega(end,2)) / 2; omega(end,end) = (omega(end-1,end) + omega(end,end-1)) / 2; end

函数里的下标关系容易看晕，所以我在关键行都写了注释。dv/dx用v在x方向中心差分，du/dy用u在y方向中心差分，两者相减，正负号千万别搞反。我早期就栽过，最终出来的图正负涡量完全互换，还盯着图看了半天没发现。

入口出口的x方向边界处理，实际要看流场数据的性质。如果数据来自周期性流向的DNS，那应该做周期边界重构；如果来自实验段中部，用零梯度外插即可。示例代码里的零梯度外插，对这种可视化目的完全够用。

3.3 主流程的调用和使用

主脚本调用方式非常简单，但有一行值得单独说：如果数据噪声大，先对速度场做一编高斯平滑，不要直接计算涡量。

% 三维网格变换后，涡量分辨率直接受速度噪声影响 u_smooth = imgaussfilt(u, 0.8); % 高斯平滑，sigma=0.8个网格 v_smooth = imgaussfilt(v, 0.8); omega = vorticity2d(u_smooth, v_smooth, dx, dy); % 如果不需要平滑，直接用原始速度 % omega = vorticity2d(u, v, dx, dy);

注意：imgaussfilt用的sigma以网格点数为单位，不是物理单位。对127x127的PIV网格，sigma在0.5到1.5之间比较合适，太小没用，太大会把细碎的涡结构全部抹平。这个参数是“平滑强度”与“结构保持”之间的平衡杆，每次处理真实数据都应该做一组敏感性测试，而不是固定照搬。

3.4 瞬时涡量场的可视化出图

算法算完了，重头戏在出图。瞬时涡量场常用的呈现方式有以下几种：

1. 涡量云图（contourf / pcolor）：底色填充，直观反映涡量空间分布
2. 等值线（contour）：叠加等值线，突出梯度变化
3. 速度矢量（quiver）：叠加在云图上，展示流动方向
4. 涡量剖面线（plot）：抽取某一x位置或y位置的涡量曲线

figure('Color', 'w', 'Position', [100 100 900 360]); % 对称颜色轴，正负涡量用红蓝区分 cmap = bluewhitered(64); % 需要自带函数，或者用flipud(jet)加白色 colormap(cmap); imagesc(x, y, omega); axis xy; axis equal; axis tight; hold on; % 叠加速度矢量，稀疏显示避免太密 skip = 4; quiver(x(1:skip:end), y(1:skip:end), ... u(1:skip:end, 1:skip:end), v(1:skip:end, 1:skip:end), ... 1.2, 'k', 'LineWidth', 0.4); caxis([-max(abs(omega(:))) max(abs(omega(:)))]); colorbar; xlabel('x / H'); ylabel('y / H'); title(sprintf('Instantaneous vorticity field, t = 1.0')); set(gca, 'FontName', 'Helvetica', 'FontSize', 12);

这里关键的是caxis([-max(abs(omega(:))) max(abs(omega(:)))])这一行。涡量有正有负，如果用默认的caxis会让零涡量不落在色轴中央，图上同一个颜色代表的意义就会发生偏移，看着像全场都是正涡量，非常误导。强制对称色轴之后，红色代表正向旋转，蓝色代表负向旋转，白色或者浅色代表零涡量区，一眼就能判断旋转方向。

如果用的是R2022a之后的MATLAB版本，官方推荐用clim替代caxis，但为了兼容旧脚本，caxis目前仍可用。建议自己写一行兼容代码：

if exist('clim', 'file') lim = [-max(abs(omega(:))) max(abs(omega(:)))]; clim(lim); else caxis([-max(abs(omega(:))) max(abs(omega(:)))]); end

这样新旧版本都不会报错。

4. 可视化参数与常见坑位

4.1 用对称色轴让正负涡量不打架

上一节已经写了对称色轴的核心代码，这里把原理讲透。涡量场包含大量零值区域，如果色轴中心不是零，深红色可能表示一个中等正涡量，深蓝色也可能表示一个更小的正涡量——颜色不再代表“旋转方向”，而是代表“相对偏离中心值的程度”，图就废了。

对称色轴还有个附加好处：在涡量云图上叠加等值线时，零涡量等值线会自动成为分界线，视觉上把正向涡核和负向涡核分开。这是排查数据质量的好帮手——如果零线穿过了原本应该平滑的剪切层内部，要么差分格式有问题，要么数据噪声没有被平滑充分。

4.2 等值线级数的选择别迷信默认值

MATLAB的contour默认自动选择约10级等值线，但20个色带在瞬时涡量场里往往会掩盖细节。我自己比较偏好的做法是先画出全幅涡量云图，观察量级分布，再手动指定等值线级数。

levels = linspace(-max(abs(omega(:))), max(abs(omega(:))), 21); contour(x, y, omega, levels, 'k', 'LineWidth', 0.3);

21级等值线，奇数保证了0等值线正中间一个。等值线太密（比如超过50级）会让云图变成一片黑网，太疏（少于9级）又看不到涡核边界。0.3线宽是出印刷级图的一个推荐值，能在多数期刊图片里保持清晰又不喧宾夺主。

4.3 速度矢量叠加的疏密与缩放

叠加quiver是最容易把图画乱的操作。通道流基流沿法向梯度很大，直接把每个网格的速度矢量画出来，近壁区域箭头密密麻麻，中心区域稀疏，视觉上严重失衡。上文代码里用了skip=4每4个网格画一个矢量，但更精细的做法是按涡量空间尺度来确定抽稀步长，保证一个涡核内部能看到至少3~5个速度矢量。

quiver的最后一个参数1.2是箭头自动缩放因子。默认值1可能导致箭头过长串出整幅图，0.5到1.5之间都需要手动试一张。打印成PDF后，箭头长度过小会糊成一片，最好用矢量格式输出时再微调一次。

4.4 网格数据不干净引起的伪涡量

真实PIV数据和粗糙CFD插值结果都免不了局部坏点或尖锐跳跃。这些坏点在做涡量差分时会被二次放大，直接表现为云图上孤立的极正或极负的点。处理这类问题，除非坏点比例极大（超过5%），否则我建议先做速度场的异常值剔除和线性插值，再做高斯平滑，最后计算涡量。

如果固定忽略一步，直接对含尖峰的速度场做差分，那出来的涡量云图几乎必然需要反复重画。动手处理前，多花两分钟用surf或pcolor扫一眼速度场本身有没有明显噪点，能省掉后面整套调试的力气。

5. 一个从解析到实测的推荐工作流

对刚入手的同学，我给一个确定性较高的流程，跑通后再去处理自己的数据：

1. 用3.1中的解析流场构造代码，生成瞬时通道流速度场；
2. 调用vorticity2d计算涡量；
3. 用3.4的可视化脚本出图，观察两个高斯涡的位置：流向三分之一处一个正向涡，三分之二处一个负向涡，近壁区应有基流贡献的强涡量带；
4. 对比解析涡量（把基流涡量 -2U0 y/H^2 与高斯扰动涡量解析表达式相加），计算全场的最大绝对误差，验证算法；
5. 确认算法无误后，替换成自己的实测速度数据，必要时加入imgaussfilt平滑。

在MATLAB里跑通全流程大约只需要十分钟。我第一次折腾这个场景时，卡在边界单边差分那里大半天，总想着把边界点去掉算了，结果画出来的图上下各有一条跟主流方向平行的异常色带，后来才发现是边界没处理。如果你也出现类似现象，第一条要检查的就是边界行和边界列到底参与了差分没有，而不是怀疑可视化参数。

另一个值得单独拎出来的细节是量化误差检验。当你用解析流场做验证时，数值涡量与解析涡量之差的绝对量级，应该在差分截断误差允许范围内。如果误差远大于理论预期，别急着观察云图结构，回头检查网格步长和差分公式的匹配，确定步长到底用的是“网格索引间隔”还是“物理坐标间隔”。gradient函数第二个参数是坐标步长不是网格点数，我见过不少人在这一步栽跟头，出来的涡量整体大一个量级，还以为是色轴标错了。