简介这套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, FPOINT\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 末尾的 FPOINT 表示数据按逐点方式排列而非块状排列选错会导致 Tecplot 把整块数据读成乱序。这个文件如果网格规模在 200 万点以上ASCII 格式会膨胀到几百 MB读取很慢。这时候可以改用 Tecplot 的二进制格式或者先把数据写成分块排布的 FBLOCK字段按维度排体积能小不少。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 里设置方法是在等值面面板中添加多个等值面值然后分别调整透明度。这套显示参数组合适合大多数 FPOINT 的稳态场数据。瞬态场数据则不建议把每一帧都导成 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涡核半径 rc5mm最大切向速度 Uθmax1m/s涡轴沿 z 方向。把解析速度场导入脚本输出 Liutex 向量后检查两个硬指标。第一个指标是旋转轴方向理论解里涡轴就是 z 方向所以 Liutex 向量的 x、y 分量应该接近零z 分量应该接近正值。在涡核中心那一列点上打印出 LiutexZ如果符号一致、数量级和理论角速度对得上说明特征分解和方向校正部分正常。第二个指标是旋转强度峰值R 的理论峰值在涡核半径处大约为 2 × Uθmax / rc 400 s⁻¹实际算出来在 380~420 之间都算正常超出这个范围就得回头查虚部阈值和梯度差分步长。验证通过后再跑一个三维倾斜涡的案例把涡轴旋转 45 度检查 Liutex 向量的三个分量是否跟随涡轴一起偏转。这一步能排查旋转轴方向标定的问题确认无误后就可以放心处理真实流场。最后给你一个我踩过坑之后养成的习惯每处理一个新案例我不会直接出图而是先把 R 的最大值、涡核中心点坐标和趋向性打印出来核对速度量级和网格尺寸是否匹配预期。这些数据不会骗人比看图猜结构可靠得多。希望这次关于 Liutex 的拆解能帮你少走些弯路拿到资源后务必备好解析验证的数据集再动手进入自己的项目。本文还有配套的精品资源点击获取