简介：这套Liutex（第三代涡识别方法）的Matlab-Tecplot实现资料，面向流体力学、CFD后处理及涡识别方向的研究者和工程师，帮助理解Liutex理论并快速上手代码实现。资源共6个文件，以3个Matlab脚本（负责核心计算、矩阵处理与Liutex导出）为主，另含PDF文献、PPT方法介绍及eqn公式文件，压缩包仅6.67MB，小巧但知识点密集。内容从涡定义讲到实际计算流程，读者可对照文档理解Liutex的张量运算，再利用现成脚本完成数据输出与Tecplot可视化，节省从零写代码和调试的时间。附带的PDF与PPT梳理了Liutex的数学背景和涡识别发展脉络，eqn公式文件便于对照关键推导，使这套包不仅是代码工具，也是很好的理论入门材料。目前已有3582人学习，适合需要系统掌握第三代涡识别方法、做科研复现或课程教学的高年级本科生、研究生与科研人员；资源免费开放，无需积分即可下载，进一步降低了学习门槛。
1. Liutex 不是又一种涡判据：先搞清楚它解决什么问题
做流体后处理时最头疼的，从来都不是算得慢，而是算出来的东西说不清。用涡量显示旋涡，剪切层里全是假阳性；用 Q 判据画等值面，阈值要靠肉眼调半天，而且不同工况之间没法横向比。这些传统判据都卡在同一个问题上：它们描述的“涡”不是直观意义上的旋转，而是数学上的某个张量特性。Liutex 方法把这个问题换了个问法——不再问“这里像不像有涡”，而是直接算“这个点的局部旋转轴指向哪个方向、旋转强度是多少”，输出的是一个三维向量场，后续无论是切片的涡核线还是三维等值面，物理意义都更明确。这份资源给的是 Liutex 理论在 Matlab 端从速度场到 Liutex 向量的完整实现，以及把结果导入 Tecplot 做可视化的配套脚本，适合正在写论文需要一种可靠涡识别方法、或者做 CFX/开源求解器后处理想换个判据验证结论的人。
2. Liutex 的数学骨架：旋转轴、旋转强度与特征分解
2.1 为什么传统判据说不清“旋转轴”方向
先对齐一个基本概念：判断一个局部流动是否在旋转，不能只看速度本身，而是看速度的空间变化率。这个变化率由速度梯度张量描述，九个分量排成一个 3×3 矩阵：
∇V = [du/dx du/dy du/dz] [dv/dx dv/dy dv/dz] [dw/dx dw/dy dw/dz]它可以分解成对称部分（应变率张量 S）和反对称部分（旋转张量 Ω）。Q 判据看的是 ‖Ω‖² - ‖S‖² 的正负，λ2 判据看的是 S² + Ω² 的特征值，这两种方法本质上都在做同一件事：拿“剪切/拉伸”和“旋转”比谁更占优，而不是把旋转单独提取出来。所以同样的流动，阈值取 0.5 和取 1 画出来的拓扑结构完全不一样，这就是玄学来源。
Liutex 的关键变化是放弃这种比例比较，直接对 ∇V 做特征分解，找出张量中代表真实旋转的分量。数学上有三种情况：一种是三个实特征值，代表该点没有局部旋转（纯变形）；另一种是一个实特征值加一对共轭复特征值，这对复特征值的存在就意味着该点存在局部旋转。Liutex 认为，复特征值的虚部决定了旋转的强度，而实特征值对应的特征向量方向，就是局部旋转轴的方向。这个思路比在一堆张量范数里找阈值更贴近物理直觉。
2.2 从特征值到旋转强度的标定关系
具体算的时候，会得到如下关系。设 ∇V 的特征值为 λr 和 λcr ± iλci，其中 λr 是实数，λcr 和 λci 是共轭复特征值的实部和虚部。Liutex 向量定义为：
R = R * r其中 r 是实特征值 λr 对应的单位特征向量（也就是旋转轴方向），而旋转强度 R 的大小取 2|λci|。为什么是 2 倍？因为在纯旋转流动里，速度梯度张量的复特征值虚部正好等于角速度的一半，乘 2 之后，Liutex 强度的量纲才能和涡量 ω 对齐（单位都是 1/s），这样方便和后处理里的其他物理量横向比较。
这里有一个值得注意的细节：当三个特征值全是实数时，表示该点没有局部旋转，直接把 Liutex 设为零向量。实际计算时，还会遇到数值误差导致虚部为微小值的伪复特征根，需要设置一个虚部阈值，小于该阈值的按实特征值处理。我在实现时用的阈值是 1e-10，但这取决于数据精度，如果是单精度浮点场可能需要放大到 1e-6 级别。
对比起来，Q 判据需要取阈值，λ2 判据也需要取阈值，而 Liutex 的旋转强度是逐点计算出来的，不需要人为设阈值。这并不代表 Liutex 在可视化阶段不需要选等值面，但至少“旋转强度”本身是客观的，不是通过一个方程左右两边谁大的不等式判出来的。这也是这个资源里最值钱的部分：它把计算核心从“猜阈值”变成了“解特征值”。
3. Matlab 端实现：从速度场到 Liutex 向量的完整流程
3.1 数据准备与速度梯度计算
这套实现接收的是结构化网格上的速度场，三个速度分量分别存成三维数组 U、V、W，维度顺序建议是 (nx, ny, nz)，对应物理坐标 (x, y, z)。如果手里是 CFD 软件导出的非结构化网格，需要先插值到矩形网格上再喂给脚本，这一点后面避坑章还会专门说。
速度梯度的计算是整套流程的精度基础。内部点用中心差分，边界用单侧差分。我这里给出一个可用的实现片段：
function [dudx, dudy, dudz] = calc_gradient(U, dx, dy, dz) [nx, ny, nz] = size(U); dudx = zeros(nx, ny, nz); dudy = zeros(nx, ny, nz); dudz = zeros(nx, ny, nz); % 内部中心差分 dudx(2:end-1, :, :) = (U(3:end, :, :) - U(1:end-2, :, :)) / (2*dx); dudy(:, 2:end-1, :) = (U(:, 3:end, :) - U(:, 1:end-2, :)) / (2*dy); dudz(:, :, 2:end-1) = (U(:, :, 3:end) - U(:, :, 1:end-2)) / (2*dz); % x 方向边界：一阶单侧差分 dudx(1, :, :) = (-3*U(1,:,:) + 4*U(2,:,:) - U(3,:,:)) / (2*dx); dudx(end, :, :) = (3*U(end,:,:) - 4*U(end-1,:,:) + U(end-2,:,:)) / (2*dx); % y、z 边界同理，此处省略 enddy 和 dz 方向的边界处理与 x 方向完全一致，不再重复展开。中心差分的截断误差是二阶，对大多数仿真数据够用；如果网格长宽比很大（比如边界层里法向网格极密），建议对法向梯度改用加权差分，否则切向梯度会被噪声放大。
这个函数返回的是 U 分量三个方向的偏导。为了组装完整的梯度张量，需要分别对 V 和 W 分量也调用一次，然后把九个偏导按矩阵排列。检查梯度的正确性有一个土办法：取一个已知旋转速度场（后面验证章会给），算出来的涡量 ω = ∇×V 应该和分析值吻合，误差在 1% 以内再继续，不然后面全白做。
3.2 逐点特征分解与 Liutex 向量组装
得到九个梯度分量之后，在每个网格点上组装梯度张量，然后做特征分解。这里我给出核心循环的完整逻辑：
% 预分配输出数组 LiutexX = zeros(nx, ny, nz); LiutexY = zeros(nx, ny, nz); LiutexZ = zeros(nx, ny, nz); omega_x = dWdy - dVdz; % 涡量分量，用于方向矫正 omega_y = dUdz - dWdx; omega_z = dVdx - dUdy; for k = 2:nz-1 for j = 2:ny-1 for i = 2:nx-1 % 组装速度梯度张量 GradV = [ dUdX(i,j,k), dUdY(i,j,k), dUdZ(i,j,k); dVdX(i,j,k), dVdY(i,j,k), dVdZ(i,j,k); dWdX(i,j,k), dWdY(i,j,k), dWdZ(i,j,k) ]; % 特征值分解 [eigVec, eigVal] = eig(GradV); lambda = diag(eigVal); % 找实特征值：虚部绝对值小于阈值 isReal = abs(imag(lambda)) < 1e-10; if sum(isReal) == 1 realIdx = find(isReal, 1); r = real(eigVec(:, realIdx)); r = r / norm(r); % 归一化旋转轴方向 % 找到一对共轭复特征值，取虚部 complexIdx = find(~isReal); R_mag = 2 * abs(imag(lambda(complexIdx(1)))); % 旋转轴方向与涡量方向对齐，保证右手定则 omega = [omega_x(i,j,k), omega_y(i,j,k), omega_z(i,j,k)]; if dot(r, omega) < 0 r = -r; end LiutexX(i,j,k) = R_mag * r(1); LiutexY(i,j,k) = R_mag * r(2); LiutexZ(i,j,k) = R_mag * r(3); end % 三个实特征值：无旋转，保持零向量 end end end特征分解用的是 Matlab 自带的 eig 函数，对小规模网格速度足够快，但要注意：如果网格是千万级体量，逐点用 eig 在循环里跑会非常慢，一个常见的优化是改用 LAPACK 的 dgeev 接口或者把张量批量展开后向量化，速度能提一个数量级。资源里给的是逐点版本，优先保证清晰度，性能不够时再改。
旋转轴方向与涡量方向对齐这一步非常关键。eig 返回的特征向量方向在数值上是任意的（正负都可能），同一物理点在不同满足度下可能得到完全相反的旋转轴方向。用涡量方向做矫正，能确保相邻网格点上旋转轴方向连续变化，画出来不会出现轴方向突变导致的涡核线断裂。
关于旋转强度的符号：Liutex 方法中 R 本身定义为非负，方向由旋转轴向量承担，所以这里取了虚部的绝对值。如果你习惯标量显示旋转强度，可以只取 R_mag，而向量方向用 r 另外存储。这两种用法在资源里分别对应了 Tecplot 的“云图显示”和“矢量显示”两种模式。
3.3 边界处理与相邻点连续性检查
上面的循环只处理了内部点，边界上的 Liutex 值沿用零或者用单侧差分补算，这一点要看实际需求。如果涡核不靠近边界，直接设零不影响显示；但如果关注的是管道近壁区域，就必须在边界上用单侧梯度补算 Liutex。补算时需要注意，边界上特征分解结果往往很不稳定，因为单侧差分精度只有一阶，张量的实特征值会频繁出现微小虚部，虚部阈值在这里要适当放松到 1e-6 级别。
做完逐点计算，我建议顺手检查一次旋转轴方向的连续性。方法很简单：取涡核中心线附近的一排点，把 Liutex 向量的三个分量打印出来，观察其方向是否呈平滑变化。如果出现间隔点方向翻转，多半是涡量方向矫正在零涡量点附近失效，可以再加一个约束：当 |ω| 小于某个阈值（比如最大涡量的 1%）时，沿用上一个网格点的旋转轴方向；这个修正能消除大部分涡核线上的毛刺。
4. Tecplot 可视化：把 Liutex 结果变成能写进论文的图
4.1 从 Matlab 导出 Tecplot 可读的 ASCII 数据
Matlab 算出的 Liutex 三个分量要和原始速度场一起导出，方便在 Tecplot 里做多变量对照。最省事的方式是写一个标准的 Tecplot ASCII 文件，点格式，坐标在前，速度次之，Liutex 分量放在最后。代码块如下：
filename = 'liutex_export.dat'; fid = fopen(filename, 'w'); fprintf(fid, 'TITLE = "Liutex Field"\n'); fprintf(fid, 'VARIABLES = "X", "Y", "Z", "U", "V", "W", "LiutexX", "LiutexY", "LiutexZ"\n'); fprintf(fid, 'ZONE I=%d, J=%d, K=%d, F=POINT\n', nx, ny, nz); for k = 1:nz for j = 1:ny for i = 1:nx fprintf(fid, '%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n', ... X(i,j,k), Y(i,j,k), Z(i,j,k), ... U(i,j,k), V(i,j,k), W(i,j,k), ... LiutexX(i,j,k), LiutexY(i,j,k), LiutexZ(i,j,k)); end end end fclose(fid);文件头里最关键的两行是 VARIABLES 和 ZONE。VARIABLES 的变量名必须和后面数据列一一对应，Tecplot 加载时按字符串匹配变量名；ZONE 末尾的 F=POINT 表示数据按逐点方式排列而非块状排列，选错会导致 Tecplot 把整块数据读成乱序。这个文件如果网格规模在 200 万点以上，ASCII 格式会膨胀到几百 MB，读取很慢。这时候可以改用 Tecplot 的二进制格式，或者先把数据写成分块排布的 F=BLOCK，字段按维度排，体积能小不少。
4.2 云图、等值面与切片：参数怎么定
数据加载进 Tecplot 后，左侧 Zone 面板里能直接看到新增的变量。显示 Liutex 强度用等值面最直观，因为旋转强度的数值是客观算出来的，不需要像 Q 判据那样纠结阈值。具体取值我一般这样定：先用菜单 Compute → 计算 Liutex 向量的模长，放到一个新变量 R 里，点击 Histogram 看分布，再选一个对应涡核区域占总体积大约 5%~10% 的值作为等值面阈值。对于典型的圆柱绕流或者机翼尾迹，R 的值通常在 2~20 之间，低于 1 的区域基本都是剪切噪声。
插一句切片问题：直接在流场中间切一个平面，用云图显示 R 的分布，颜色表选用 Tecplot 默认的“Blue-White-Red”，把最小值和最大值对标到涡核内外，这样平面图上能看到涡核的圆斑，和 Q 判据云图相比，边界更锐利。如果你需要显示旋转轴方向，把 LiutexX/Y/Z 作为矢量变量叠加到切面上，箭头方向就是局部旋转轴，注意 Tecplot 矢量图的长度比例因子要调到能区分相邻点方向的程度，通常取网格间距的 0.5 倍。
4.3 三维显示里旋转轴的呈现技巧
三维等值面只能看强度分布，看不出旋转轴的连续变化。一个有效做法是画一条过涡核中心的线，沿线采样 Liutex 方向，用 Tecplot 的 Polyline 功能从涡核入口追到出口，然后把沿线 Liutex 向量的三个分量画成曲线图。这样可以直观看到旋转轴在空间里的偏转角度，这在论文里比一张彩色等值面图更有说服力。
画三维半透明等值面时，我习惯把 R 的最小等值面阈值设为峰值 R 的 30% 到 50%，让透明层显示外圈、不透明层显示内核。在 Tecplot 里设置方法是在等值面面板中添加多个等值面值，然后分别调整透明度。这套显示参数组合适合大多数 F=POINT 的稳态场数据。瞬态场数据则不建议把每一帧都导成 ASCII 文件，最好在 Matlab 里把多帧的 Liutex 向量统一存成一个 NETCDF 文件，利用 Tecplot 对 NETCDF 的原生支持做时间动画，否则帧数一多文件管理会乱成一团。
5. Liutex 使用避坑：五个容易翻车的操作与排查
5.1 旋转轴方向频繁反转
现象：画涡核线矢量图时，相邻网格点的旋转轴方向出现 180 度翻转，显示效果像是箭头的分布随机跳变，涡核中心线上表现得尤其明显。
原因：eig 返回的特征向量方向本来就没有固定朝向，负方向也是合法解；只用涡量方向做校正，在涡核中心附近涡量很小甚至为零，dot(r, ω) 的符号被噪声主导，导致方向随机翻转。
解决：在方向校准前先用最小涡量阈值过滤，|ω| 小于最大涡量 1% 的点不翻转，直接用相邻已确定的梯度方向，让涡核内部的轴方向顺着流向上游的一致性趋势走。我实际做法是分两遍：第一遍正常算 Liutex 向量，第二遍把该层每个方向的特征向量分量的第一个非零分量符号对齐到前一个点，两遍过后轴方向基本连续。
5.2 虚部阈值设置不当引出假旋转
现象：在无旋的均匀剪切流动区域，计算出来的 Liutex 强度不为零，有时候在边界层里还能画出大片的“伪涡”。
原因：单精度速度场的梯度噪声达到 1e-6 量级，梯度张量的特征值虚部也会出现微小的伪复数值，如果阈值设成 1e-10，这些伪复特征根就会全被当成真实旋转。
解决：检查数据精度。单精度场把虚部阈值提高到 1e-6 甚至 1e-5，双精度场再逐步降低；同时结合涡量做辅助筛选，如果 Liutex 强度很大但该点涡量几乎为零，这几乎一定是数值噪声，直接剔除。
5.3 网格维度顺序导致 Tecplot 数据错位
现象：Tecplot 里云图显示时，强度分布明明是对的，但坐标轴的方向乱了，切片也切不到预期的位置。
原因：Matlab 里数组维度顺序是 (nx, ny, nz) 对应 (x, y, z)，写导出时用的是三重循环，循环顺序或索引顺序与 ZONE 里声明的 IJK 顺序不匹配，数据列当然对不上坐标。
解决：导出前把 ZONE 行里的 I、J、K 和循环外层对应好。我通常统一用最外层 k 循环对应 Z，中间 j 循环对应 Y，最内层 i 循环对应 X，强制约束住就不会错。导出后取样首尾若干行核对 X 坐标是否从最小渐变到最大，确认后再进 Tecplot。
5.4 等值面阈值选错看不出涡核结构
现象：画 R 的等值面时，阈值调来调去，要么看不到东西，要么整个流动区域全被包成一个大团块。
原因：R 是客观值但没有“归一化”，不同工况量级差很大。低速流场 R 才 0.1，高速流场 R 超过 100，拿上一步用的 5 去套新数据当然失灵。
解决：每换一个案例先算 R 的最大值和直方图，把最大值的 1/3 作为第一个等值面试。另外不要直接用 Liutex 向量分量画等值面，那些分量有正负，会同时出现正反两组面，只会把图弄乱。
5.5 非结构化网格上的梯度计算偏差
现象：同一套物理流动，从非结构网格插值到结构化网格后，Liutex 强度与原始求解器后处理结果差异很大，涡核线形状都对不上。
原因：插值过程本身会抹掉一部分速度梯度信息，特别是边界层和剪切层位置，梯度大、网格密，插值误差容易被放大。
解决：先检查原始网格在目标区域的网格密度是否和结构化网格匹配。常用做法是先把原始网格上速度场用高阶插值（比如 RBF 插值）转出来，再用最小二乘法在目标网格上重构梯度，这样比简单线性插值保真度高很多。如果两个网格差异实在太大，就别强求 Liutex 数值完全一致，重点比较涡核位置而不是强度绝对值。
6. 用最小涡验证 Liutex 实现：一个半小时能做完的自检流程
拿到这套资源后，我强烈建议你先跑一遍内置的验证案例，而不是直接拿来算自己的流场。验证流场用解析的 Lamb-Oseen 涡，二维足够，因此切成一个很小的三维网格。计算域设 60×60×3 个点，物理尺寸 0.1m × 0.1m × 0.005m，涡核半径 rc=5mm，最大切向速度 Uθmax=1m/s，涡轴沿 z 方向。把解析速度场导入脚本，输出 Liutex 向量后检查两个硬指标。
第一个指标是旋转轴方向：理论解里涡轴就是 z 方向，所以 Liutex 向量的 x、y 分量应该接近零，z 分量应该接近正值。在涡核中心那一列点上打印出 LiutexZ，如果符号一致、数量级和理论角速度对得上，说明特征分解和方向校正部分正常。第二个指标是旋转强度峰值：R 的理论峰值在涡核半径处，大约为 2 × Uθmax / rc = 400 s⁻¹，实际算出来在 380~420 之间都算正常，超出这个范围就得回头查虚部阈值和梯度差分步长。
验证通过后，再跑一个三维倾斜涡的案例，把涡轴旋转 45 度，检查 Liutex 向量的三个分量是否跟随涡轴一起偏转。这一步能排查旋转轴方向标定的问题，确认无误后就可以放心处理真实流场。
最后给你一个我踩过坑之后养成的习惯：每处理一个新案例，我不会直接出图，而是先把 R 的最大值、涡核中心点坐标和趋向性打印出来，核对速度量级和网格尺寸是否匹配预期。这些数据不会骗人，比看图猜结构可靠得多。希望这次关于 Liutex 的拆解能帮你少走些弯路，拿到资源后务必备好解析验证的数据集，再动手进入自己的项目。
本文还有配套的精品资源，点击获取
