Matlab与Tecplot实现Liutex涡识别全流程详解

简介：该资源为Liutex（涡定义与第三代涡识别方法）的入门与实现资料包，面向流体力学、涡动力学方向的研究人员及工程技术人员，帮助读者理解Liutex理论并快速上手Matlab与Tecplot的可视化实现。压缩包共6个文件，包括3个Matlab脚本（如getR.m、writeLiutex.m等，用于计算与输出Liutex场）、1份PDF文献介绍、1份PPTX讲解材料及1个公式文件，整体大小6.67MB，内容精简但覆盖从原理到编码的关键环节。已有3582人浏览学习，较为适合具有一定流体力学基础、希望掌握新代涡识别方法的读者。资料由作者针对网上盗卖行为重新免费上传，无需积分即可下载。通过阅读PDF与PPT可系统了解Liutex定义与优势，借助Matlab脚本可对照实现核心算法，并进一步结合Tecplot完成结果可视化，从而在科研或工程项目中直接应用这一较前沿的涡识别工具。

1. Liutex 不是玄学：为什么我用 Matlab 算、用 Tecplot 看

做 CFD 后处理的人大多会遇到同一种尴尬：算完流场，用 Q 准则或 λ2 准则挑涡结构，结果射流剪切层、壁面边界层、甚至纯剪切缝隙里都飘着一大片“涡”，可这些区域根本没有旋转核心，云图出来后完全没法跟实验流动显示对上。Liutex 方法就是来补这个空子的，它把速度梯度张量拆成刚性旋转和变形两类贡献，只保留真正的刚性旋转强度 R，再输出矢量。这条路我用 Matlab 做逐点计算、Tecplot 做后处理出图，脚本全流程跑通，不依赖商业 CFD 求解器内部函数。适合手上只有速度场结果、想快速验证涡结构的仿真工程师和研究生。

2. Liutex 的数学原理与计算路线：从速度梯度到旋转轴

2.1 为什么不再依赖 Q 准则和 λ2 判据

传统判据的通病，在于它们都建立在“涡量大于变形就等于有涡”这个近似上。考察速度梯度张量∇u的对称部分 S 和反对称部分 Ω，Q 准则计算0.5*(||Ω||² - ||S||²)，λ2 准则看的是 S²+Ω² 的第二特征值。这两个判据在纯剪切层里的响应非常糟糕：纯剪切场的涡量和变形量可以同时很大，Q 值照样偏高，于是把一大片剪切区误报成涡。我最早用 Q 准则时，判断“哪里是涡”基本靠猜，后来才意识到这是剪切污染。

Liutex 的计算思路完全不同。它不比较涡量与变形的相对大小，而是先问一个问题：这个局部区域是否存在一条旋转轴，能让流体质点绕它做刚性旋转？判断依据是速度梯度张量的特征值结构。只要存在一对共轭复特征值和一条唯一的实特征值对应的旋转轴，再满足强旋转判别式，就认为该点是涡。这个门槛把剪切、拉伸这类非旋转变形滤掉了，所以它输出的等值面通常比 Q 准则小两三圈，正好对应流体绕某条轴转动的核心。

2.2 特征值分解与 Liutex 强度的计算

完整算法分五步。第一步，计算每一点的速度梯度张量∇u；第二步，求∇u的三个特征值，判断是否存在一个实特征值和一对共轭复特征值，如果三个特征值都是实数，该点直接置零；第三步，取实特征值对应的特征向量，单位化后作为旋转轴方向 r；第四步，计算涡量 ω 在 r 上的投影ω_r = ω·r，并取复特征值的实部 λ_cr；第五步，求刚性旋转强度 R。公式是：

R = sign(ω_r) * sqrt( max(ω_r² - 4*λ_cr², 0) )

Liutex 向量就是R·r。公式里的4*λ_cr²项，就是把速度梯度张量中由拉伸/剪切带来的变形贡献扣除，剩下才是绕 r 轴的刚性旋转分量。这套计算在 Matlab 里可以做成一个独立函数，方便后续逐点调用。

function [Rvec, Rmag] = liutex_grad(D, omega) % D : 3x3 速度梯度张量 % omega : 1x3 涡量向量 % 返回 Liutex 向量 Rvec 与标量强度 Rmag tol = 1e-8 * max(1, sqrt(sum(D(:).^2))); [V, Deig] = eig(D); lam = diag(Deig); n_real = sum(abs(imag(lam)) < tol); Rvec = [0 0 0]; Rmag = 0; if n_real ~= 1 return; % 三个实根或退化情况，不产生刚性旋转 end idx_real = find(abs(imag(lam)) < tol, 1); r = real(V(:, idx_real)); r = r / norm(r); if dot(r, omega) < 0 r = -r; % 统一旋转轴方向，避免相邻点符号跳变 end idx_cplx = find(abs(imag(lam)) >= tol, 1); lam_cr = real(lam(idx_cplx)); wr = dot(omega, r); disc = wr * wr - 4 * lam_cr * lam_cr; Rmag = sign(wr) * sqrt(max(disc, 0)); Rvec = Rmag * r; end

逻辑说明：每个网格点单独调用这个函数，返回零向量说明该点不是涡。n_real用容差判断而不是imag(lam) == 0，因为浮点计算里纯实数根的虚部也可能出现 1e-16 量级的残留，严格相等会把好点全部漏掉。旋转轴方向统一那一步不是可选项，特征向量只有方向没有正负，同一物理涡在相邻网格点可能一个取正一个取负，不做对齐则输出 R 向量会连续倒向。

参数说明：tol取张量 F 范数的 1e-8 倍，比固定阈值如 1e-10 更抗量纲影响，压气机流场速度量级大、梯度张量数值大，固定阈值很容易把所有特征值都判成实根；lam_cr取第一个虚部超过容差的特征值的实部，一对共轭复特征值实部相同，取哪一个结果一样。R 和 λ_cr 的单位都是 1/s，与涡量一致，所以后处理时可以直接和涡量场对比量级。

2.3 旋转轴方向统一：关键数值细节

特征向量作为旋转轴方向时，存在一个天然歧义：r和-r是同一个方向。如果不在代码里做统一，输出 R 向量会在相邻网格点之间出现随机反向。我一般用dot(r, omega) < 0来判断，让 r 与局部涡量方向的点积非负。这样旋转轴方向在空间上连续，R 向量在流场中的指向才有物理意义。

极少数退化点上，涡量接近零或 r 与 ω 接近垂直，点积接近零，符号判定会抖动。这类点本身 R 也趋于零，对涡核等值面影响可以忽略。不要在代码里把这类点特殊处理成 NaN，Tecplot 显示时会留下空洞，反而干扰后续处理。

3. Matlab 实现：从 CFD 速度场到 Liutex 矢量的完整流程

3.1 读入结构化网格速度数据

假设已经有计算出的结构化网格速度场，无论是自己写的求解器还是从商业软件导出的 Tecplot ASCII 文件，这一步都先把数据读到 Matlab 四维数组中。Tecplot ASCII 的 F=POINT 格式比较简单，数据按 I 方向最快、J 方向其次、K 方向最慢排列，每个点按顺序写 X、Y、Z、U、V、W。我习惯先扫描到 ZONE 关键字再开始读数，避免不同软件导出的表头行数不一致。

function [x,y,z,u,v,w] = read_tecplot_ascii(fname, ni, nj, nk) % 读取 Tecplot ASCII 点序数据 % 变量顺序固定为 X Y Z U V W fid = fopen(fname, 'r'); offset = 0; while ~feof(fid) line = fgetl(fid); offset = offset + 1; if contains(upper(line), 'ZONE') break; end end raw = textscan(fid, '%f'); fclose(fid); data = raw{1}; npts = ni * nj * nk; data = reshape(data(1:npts*6), 6, [])'; x = reshape(data(:,1), ni, nj, nk); y = reshape(data(:,2), ni, nj, nk); z = reshape(data(:,3), ni, nj, nk); u = reshape(data(:,4), ni, nj, nk); v = reshape(data(:,5), ni, nj, nk); w = reshape(data(:,6), ni, nj, nk); end

逻辑说明：先用fgetl逐行找含 ZONE 的行，找到后文件指针已经停在数据起点，textscan从下一行开始读取所有浮点数。reshape成 6 列是为了把每个点的分量拆开，再用reshape还原成三维数组时，第一维对应 I 方向，和 Tecplot 的 F=POINT 写入顺序一致。

参数说明：ni、nj、nk必须与 Tecplot ZONE 行里的 I、J、K 一致，否则 reshape 会报错或静默错位。读取后建议先对 x、y、z 做一次单调性检查，diff(x(:,1,1))出现变号说明数据顺序理解反了。

3.2 核心函数：逐点计算 Liutex 向量

读入速度分量后，用 Matlab 自带的gradient函数计算九个梯度分量，再组装成速度梯度张量 D，调用上一章写好的liutex_grad。注意gradient的三参数版本要求网格等间距，结构网格且 dx、dy、dz 为常数时可以直接用。

dx = x(2,1,1) - x(1,1,1); dy = y(1,2,1) - y(1,1,1); dz = z(1,1,2) - z(1,1,1); [dudx, dudy, dudz] = gradient(u, dx, dy, dz); [dvdx, dvdy, dvdz] = gradient(v, dx, dy, dz); [dwdx, dwdy, dwdz] = gradient(w, dx, dy, dz); wx = dwdy - dvdz; wy = dudz - dwdx; wz = dvdx - dudy; Rx = zeros(size(u)); Ry = Rx; Rz = Rx; Rmag = Rx; for k = 1:nk for j = 1:nj for i = 1:ni D = [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)]; omega = [wx(i,j,k), wy(i,j,k), wz(i,j,k)]; [Rvec, Rm] = liutex_grad(D, omega); Rx(i,j,k) = Rvec(1); Ry(i,j,k) = Rvec(2); Rz(i,j,k) = Rvec(3); Rmag(i,j,k) = Rm; end end end

逻辑说明：速度梯度张量 D 的三行分别是 u、v、w 的三个方向导数，组装顺序要和liutex_grad里的eig(D)对应。涡量用旋度公式直接算，符号约定是右手系，如果数据来自左手坐标系，R 会整体反向，需要在读入阶段先检查坐标系。

参数说明：gradient在边界上用单侧差分，边界层第一层网格的梯度和 R 会偏大，这是后处理时最需要注意的误差来源，后面避坑章节会展开。三重循环在五百万网格点上大概要跑几分钟到十几分钟，主要开销在eig；如果网格过亿，可以把外层循环改成parfor，但liutex_grad内部不要共享变量。

3.3 导出 Tecplot 可直接打开的数据文件

计算完成后，把速度原始场和 Liutex 场一起写入新的 Tecplot ASCII 文件。这样导入 Tecplot 后，既可以同时显示原始速度云图和 Liutex 涡核，又能在同一份数据里做截面和向量图，不用来回切换文件。

fidw = fopen('liutex_field.dat', 'w'); fprintf(fidw, 'TITLE = "Liutex field"\n'); fprintf(fidw, 'VARIABLES = "X","Y","Z","U","V","W","Rx","Ry","Rz","Rmag"\n'); fprintf(fidw, 'ZONE T="Liutex", I=%d, J=%d, K=%d, F=POINT\n', ni, nj, nk); for n = 1:numel(x) fprintf(fidw, '%.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e\n', ... x(n), y(n), z(n), u(n), v(n), w(n), Rx(n), Ry(n), Rz(n), Rmag(n)); end fclose(fidw);

逻辑说明：numel(x)线性索引顺序与 Tecplot 的 F=POINT 顺序一致，因为数据写入时 I 方向最快。变量列表里 Rx、Ry、Rz 是 Liutex 向量三个分量，Rmag 是标量强度，Tecplot 里既能用 Rmag 做等值面，也能用 Rx/Ry/Rz 画旋转轴向量。

参数说明：%.6e足够后处理绘图精度，不建议用%f，因为速度量级跨度大，固定小数位会丢失小量信息。老版本 Tecplot 打开时报 zone 类型错误，可以在 ZONE 行加ZONETYPE=ORDERED。写完文件后，用load('liutex_field.dat')读回来检查第一列是否单调递增，这是最廉价的验证方式。

4. Liutex 实现的常见坑：现象、原因与排查

4.1 旋转轴符号跳变导致 R 向量方向混乱

现象：输出的 R 向量在相邻网格点之间突然完全反向，Tecplot 里向量图呈锯齿状，等值面边界出现随机毛刺。

原因：特征值分解给出的特征向量方向不确定，r和-r是同一根轴，逐点计算时不同网格点可能落到不同符号分支上。

解决：在liutex_grad里强制dot(r, omega) < 0时取反。如果做完这步还有个别点跳变，检查这些点的涡量是否接近零，接近零的区域 R 也接近零，不影响涡核判断，不用处理。

4.2 特征值虚部判断失误把剪切错认成旋涡

现象：层流库埃特流这种纯剪切场，Rmag 仍然和真实涡核同量级，等值面密密麻麻。

原因：代码里用了imag(lam) == 0判断实根，浮点误差让纯实数根的虚部出现 1e-16 量级残留，被误判成复数特征值，于是走进了旋转分支。

解决：用相对容差tol = 1e-8 * norm(D, 'fro')，并且只要n_real ~= 1就返回零向量。更保守的做法是同时要求disc > 0才给非零 R，根号内非正的区域一律按非涡处理。

4.3 sqrt 内判别式出现负值

现象：个别网格点 Rmag 突然出现 NaN，Tecplot 等值面在这些点附近断开。

原因：ω_r² - 4*λ_cr²理论上是非负的，但速度场经过插值或数据格式转换后，局部梯度会引入噪声，导致根号内出现 -1e-10 量级的负数；另一种可能是该点本身是强变形主导，两个共轭复特征值的实部过大，不满足强旋转条件。

解决：用max(disc, 0)保护，同时在脚本里统计disc < 0的点数。如果负值点占比小于总点数的 0.1%，可以直接忽略；如果超过 1%，说明输入速度场没有经过良好插值，或者网格变换时梯度算错了，不要靠改公式掩盖。

4.4 Tecplot 打开数据错位：变量数与 F=POINT 不匹配

现象：导入 liutex_field.dat 后，坐标云图出现条纹，或者某个变量显示出来的数值明显不对。

原因：fprintf 的占位符数量与 VARIABLES 声明的变量数不一致，常见的是少写一个分量或多写一个分量，Tecplot 会按固定列宽读取，错位后整个场都乱掉。

解决：检查 VARIABLES 字符串里有几个变量名，fprintf 里每行就写几个%e。我一般写完后立刻用textscan读回来，对比 Rmag 总量是否与计算时一致。另一个隐蔽问题是文件被 Excel 打开后另存，换行符和空格都会被改，Tecplot 可能读不出正确的点数，所以别让编辑器碰 .dat 文件。

4.5 边界差分误差让 R 值虚高

现象：壁面附近出现一层条带状伪涡，Rmag 在贴近边界的第一层网格上突然增大，向内却迅速下降。

原因：gradient在边界上只能用单侧差分，一阶精度比内部中心差分低一到两个数量级，速度梯度过大导致特征值误判，R 在边界被显著放大。

解决：先分清是物理边界还是人为截断边界。物理壁面处建议把计算域收缩一个网格，跳过最外层两点再统计；数据量足够时，也可以在计算前对速度场沿壁面法向做一次五点平滑，但要注意不要过度平滑导致真实涡核消失。边界层的 R 虚高是所有特征值类涡判据的共性问题，Q 准则同样存在，不是 Liutex 特有的。

5. 验证脚本正确性的三个解析场与 Tecplot 出图技巧

5.1 用解析速度场做量级验证

拿到新数据直接出图之前，花十分钟用三个解析速度场把脚本跑一遍，能避免大量后期返工。表格是三个验证用例：

刚体旋转用来验证主量级是否正确，纯剪切用来验证是否会被误判成涡，旋转加拉伸用来验证变形贡献是否被正确剥离。第三组要求 |ε| < |Ω|，否则没有刚性旋转主导区，R 应当返回零。

nx = 50; ny = 50; x = linspace(-1, 1, nx); y = linspace(-1, 1, ny); [X,Y] = meshgrid(x, y); Z = zeros(size(X)); Omega = 1.0; Epsilon = 0.5; u = -Omega*Y + Epsilon*X; v = Omega*X - Epsilon*Y; w = 0*Z; % 后续复用 3.2 节的梯度与循环代码 % 取中心点 (nx/2, ny/2) 的 Rmag 与 2*sqrt(Omega^2-Epsilon^2) 对比

逻辑说明：网格中心点远离边界，差分误差最小，适合做定量对比；期望值与数值解误差应小于 1e-8。如果误差很大，优先检查组装 D 时梯度分量放的行列位置是否反了，这个错误会把旋转方向算错但不影响特征值实部，非常隐蔽。

参数说明：Ω 取 1.0 是为了让量级便于检查，刚体旋转的涡量是 2Ω，所以 Rmag 期望值是 2.0，不需要额外换算。

5.2 Tecplot 里画涡核等值面和旋转轴向量

数据导入 Tecplot 后，我一般先做一个 Rmag 的等值面，阈值取全场最大值的一半，快速看涡核空间位置；再用 Rx、Ry、Rz 画 Vector，确认旋转轴方向与涡量方向大致一致。等值面不要直接拖 Contour 默认值，先看 Rmag 分布直方图，找到有明显拐点的阈值，避免被边界伪涡拉高最大值后等值面什么都显示不出来。

一个实用的步骤是：在 Tecplot 里创建两个切片，一个垂直于旋转轴看截面涡核，一个沿旋转轴看延展长度。如果涡核在切片上呈对称圆形，说明旋转轴估计正确；如果被拉成椭圆形，说明该区域存在较强拉伸，需要回到 Liutex 向量而不是只盯着 Rmag。

有一次我图省事，拿到新的模拟数据直接跑脚本出等值面，结果壁面全是伪涡结构。排查了半天才发现是特征值判断没加容差，把剪切区的实数根误判成了复共轭根。从那以后我每次换网格格式、换求解器，都强制先用这三个解析场把脚本过一遍，再上真实数据。希望帮到你。

本文还有配套的精品资源，点击获取