简介：一套围绕Liutex涡识别方法的资料包，系统提供方法介绍与对应的Matlab程序实现，适合流体力学方向的研究生、CFD工程技术人员以及对第三代涡识别方法感兴趣的科研人员。Liutex作为近年提出的涡定义，旨在克服传统涡判据在剪切流中的误判问题，这套资料恰好覆盖了从理论基础到代码落地的完整链路。包内共6个文件，包括3个.m源程序（getR.m、writeLiutex.m、partrow.m），覆盖Liutex量计算、行数据处理与后处理调用等环节；1个.pptx演示文稿和1个PDF文档对Liutex涡定义原理及实现方法进行系统讲解；另附1个.eqn公式文件，便于推导和排版时引用核心公式。压缩包整体6.67MB，结构清晰、下载快捷，目前已有3533人学习使用。读者既能借助讲义理解Liutex的数学定义与涡核识别流程，也能直接运行源代码，在Matlab与Tecplot环境下完成实际流场数据的Liutex计算与可视化，节省自行查阅文献、摸索代码的时间与精力，对入门Liutex方法和开展具体分析都有直接帮助。 我记得第一次把Liutex方法用在圆柱绕流的LES数据上时，审稿人问了一个让我很窘迫的问题：“你的Q准则等值面阈值是怎么确定的？”当时我只能含糊带过，因为Q准则这类涡识别方法本质上依赖速度梯度的量级，阈值一旦变了，涡核位置和范围都跟着漂。后来我把Liutex方法引入分析流程，用Matlab算出Liutex向量场，再用Tecplot做等值面和涡核线，才真正把“刚体旋转”和“纯剪切”分离开。这篇文章就是我从算法原理、Matlab代码实现到Tecplot出图全过程的记录，适合流体力学方向的研究生和工程师，尤其是有CFD数据但不知道怎么提取涡结构的人。
1. 为什么我把涡识别从Q准则换成Liutex：刚体旋转与剪切的分离
1.1 Q准则、lambda2和Omega的不良体感
先说实话，Q准则不是不能用，而是用起来太“手滑”。它的物理含义是速度梯度张量的二阶不变量：Q = 0.5 (||Ω||^2 - ||S||^2)，其中Ω是旋转张量，S是应变率张量。涡的区域被定义为Q > 0，也就是局部旋转占主导。问题在于，Q会同时捕捉纯剪切层的旋转，比如边界层内的高剪切区域经常被误判成涡，导致等值面“糊成一片”。
lambda2方法稍有改善，它的思路是找压力极小值条件的替代项，从对称张量 S^2 + Ω^2 的第二大特征值入手。但lambda2同样对剪切敏感，而且在低速流动中阈值很难标定。Omega方法引入了一个旋转权重系数，等于涡量平方除以涡量平方加应变率平方，虽然敏感度降下来了，但添加的小量epsilon需要人为设定，本质还是“调参”。
1.2 Liutex的定义逻辑：用实特征值定量刚体旋转
Liutex的核心思想非常直接：流场中某一点的旋转，本质上是流体质点绕某个局部转轴做刚体旋转，而不是剪切形变。这个思想通过速度梯度张量的Schur分解来实现。对三维速度梯度张量∇V做Schur分解后，可以得到一个实特征值λcr和对应的单位特征向量r。这里的λcr对应刚体旋转的旋转强度正交项（数值上关联旋转角速度的量级），r就是局部旋转轴的单位方向。
于是Liutex向量定义为：
R = λcr · r
这看起来简单，但和Q、lambda2有个本质区别：Liutex输出的是一个有方向的向量，而不只是一个标量阈值场。你不仅能判断“这里有没有涡”，还能知道涡的转轴朝哪、旋转强弱多大。这个特性在提取三维涡核线时特别有用，因为可以直接沿着r方向追踪。
1.3 Liutex带来的实际收益
我用Liutex重构圆柱绕流尾涡时，最直观的感受是涡核位置比Q准则干净得多。Q准则下尾涡区域附近总有一层高Q值的剪切层晕影，Liutex等值面则是清晰的管状结构，涡核线能直接穿过等值面中心。后续做涡量输运分解时，Liutex还能把刚性旋转对涡量的贡献单独拆出来，这在研究涡拉伸、涡撕裂时是Q准则完全做不到的。
2. Matlab里算Liutex的完整套路：梯度、Schur分解与旋转轴符号修正
2.1 速度梯度张量：gradient函数是主力，但要小心方向
Matlab计算Liutex的第一个关键步骤是速度梯度张量。对三维有序网格数据，我通常直接用gradient函数：
[dUx, dUy, dUz] = gradient(U, x, y, z); [dVx, dVy, dVz] = gradient(V, x, y, z); [dWx, dWy, dWz] = gradient(W, x, y, z);这里有个非常容易踩的坑：gradient函数处理均匀网格和理想坐标向量没问题，但如果你的网格是变间距的，务必把x、y、z三个坐标向量完整传进去，否则梯度算出来整体偏大或偏小，Liutex的幅值就会失真。另一个坑是坐标系顺序，Matlab的数组维度顺序是(行, 列, 页)，也就是(y方向, x方向, z方向)，和Tecplot常见的(I, J, K)顺序不同。我建议在写脚本时第一时间把维度顺序用注释固定下来，不然后面导出Tecplot格式能绕晕。
2.2 Schur分解与实特征值提取：eig比schur更好用
理论上的Schur分解对应标准算法schur(Grad)，但实际操作中我发现直接用eig(Grad)反而更顺手。因为速度梯度张量是一个3×3矩阵，通常有三个特征值：一个实特征值λcr，一对共轭复特征值。我要做的就是找出虚部绝对值最小的那个特征值，取它的实部作为λcr，对应的特征向量作为旋转轴r。
代码可以这样写：
% Grad 是当前网格点的 3x3 速度梯度张量 [Vec, Val] = eig(Grad); lambda = diag(Val); [~, idx] = min(abs(imag(lambda))); lambda_cr = real(lambda(idx)); r = real(Vec(:, idx)); r = r / norm(r); % 单位化注意，如果流场中存在退化点，三个特征值可能全部接近实数，这时候取虚部最小会导致误判。我通常额外加一个判断：如果计算得到的abs(lambda_cr)小于全场统计值的一个极小量（比如最大值的1e-6），就把该点Liutex设为0，避免噪声放大。
2.3 旋转轴符号歧义：必须做连续性修正
这一步是很多初稿代码里最容易缺失的部分，也是我认为Liutex实现里最需要经验支撑的地方。特征向量r在数学定义上有符号歧义，因为如果r是特征向量，那-r同样满足特征方程。如果每个网格点独立选特征向量，相邻点的r方向可能突然反转，等值面和涡核线会呈现锯齿状。
我的解决思路是连续性约束：从流场中心或边界的一个种子点出发，对周围点的r做符号对齐。简单实现是遍历相邻点，如果dot(r_current, r_neighbor) < 0，就把neighbor的r取反。对大规模数据可以按z轴分层做，先固定第一层，再沿流向逐层传播。
2.4 三维Liutex场计算的完整函数
把上面的逻辑串起来，一个可直接调用的三维场计算函数如下：
function R = liutex_field(X, Y, Z, U, V, W) % 输入: X,Y,Z 网格坐标; U,V,W 三方向速度分量 % 输出: R 为 Nx x Ny x Nz x 3 的 Liutex 向量场 [dUx, dUy, dUz] = gradient(U, X(1,:,1), Y(:,1,1), Z(1,1,:)); [dVx, dVy, dVz] = gradient(V, X(1,:,1), Y(:,1,1), Z(1,1,:)); [dWx, dWy, dWz] = gradient(W, X(1,:,1), Y(:,1,1), Z(1,1,:)); [nY, nX, nZ] = size(U); R = zeros(nY, nX, nZ, 3); threshold = 1e-6 * max(abs(dUx(:))); for k = 1:nZ for j = 1:nY for i = 1:nX Grad = [dUx(j,i,k), dUy(j,i,k), dUz(j,i,k); ... dVx(j,i,k), dVy(j,i,k), dVz(j,i,k); ... dWx(j,i,k), dWy(j,i,k), dWz(j,i,k)]; [Vec, Val] = eig(Grad); lambda = diag(Val); [~, idx] = min(abs(imag(lambda))); lambda_cr = real(lambda(idx)); if abs(lambda_cr) < threshold continue; end rVec = real(Vec(:, idx)); rVec = rVec / norm(rVec); R(j,i,k,:) = lambda_cr * rVec; end end end % 符号一致性修正: 按 X 方向逐列传播 for k = 2:nZ for j = 2:nY for i = 2:nX rPrev = squeeze(R(j,i-1,k,:)); rCurr = squeeze(R(j,i,k,:)); if dot(rPrev, rCurr) < 0 R(j,i,k,:) = -R(j,i,k,:); end end end end end这段代码没做什么花哨优化，但对中等规模网格（比如500万点以内）已经能接受。再大的网格建议改成parfor或在循环里逐层处理，避免一次性把所有梯度张量存进内存导致内存爆炸。
3. 从Matlab到Tecplot：有序网格数据的导出与导入陷阱
3.1 Tecplot能直接识别的文本格式
Matlab算完Liutex向量场之后，最稳妥的方式是导出Tecplot可读的ASCII.dat文件。虽然最新版Tecplot支持NetCDF和HDF5，但ASCII格式最通用、兼容性最好，也方便中途用记事本排查数据问题。
对有序网格，一个标准的Tecplot文件结构是：
TITLE = "Liutex field" VARIABLES = "X", "Y", "Z", "U", "V", "W", "Lx", "Ly", "Lz", "LiutexMag" ZONE I=128, J=64, K=32, F=POINT x y z u v w Lx Ly Lz LiutexMag ...其中I对应X方向格点数，J对应Y方向格点数，K对应Z方向格点数。如果你的Matlab数组第一维是Y方向，那么写文件时循环嵌套必须按 K -> J -> I 的顺序输出，我一开始就是没注意这一点，导入Tecplot后等值面全是扭曲的。F=POINT表示数据按逐点顺序排列，这是最不会出错的格式；计算精度高的话用%g格式化输出，文件大小可控。
3.2 Matlab导出脚本的公开处刑时刻
下面这段是我常用的导出脚本核心片段：
function write_tecplot_dat(filename, X, Y, Z, U, V, W, R) [nY, nX, nZ] = size(U); fid = fopen(filename, 'w'); fprintf(fid, 'TITLE = "Liutex field"\n'); fprintf(fid, 'VARIABLES = "X", "Y", "Z", "U", "V", "W", "Lx", "Ly", "Lz", "LiutexMag"\n'); fprintf(fid, 'ZONE I=%d, J=%d, K=%d, F=POINT\n', nX, nY, nZ); Lmag = sqrt(R(:,:,:,1).^2 + R(:,:,:,2).^2 + R(:,:,:,3).^2); for k = 1:nZ for j = 1:nY for i = 1:nX fprintf(fid, '%g %g %g %g %g %g %g %g %g %g\n', ... X(j,i,k), Y(j,i,k), Z(j,i,k), ... U(j,i,k), V(j,i,k), W(j,i,k), ... R(j,i,k,1), R(j,i,k,2), R(j,i,k,3), ... Lmag(j,i,k)); end end end fclose(fid); end有一点要注意：ZONE行里的 I= 后面是X方向格点数，不是Matlab矩阵第二维吗？对，因为第二维是X方向，所以nX。如果你像我一样用[nY, nX, nZ] = size(U)这种读法，逻辑就很清楚。
3.3 多块或多时刻数据的处理习惯
如果CFD数据是多块网格，不要试图放在同一个ZONE里，Tecplot无法直接处理共享边界重叠的多块ASCII数据。更稳妥的办法是每个块写成一个ZONE，同时在VARIABLES里加一个辅助变量blockID，导入后在Tecplot里用blockID做条件过滤。多时刻数据则建议每个时刻单独一个.dat文件，然后用Tecplot的批量导入功能，一次加载整个序列做动画。
4. Tecplot可视化三板斧：等值面、涡核线、高清图导出参数
4.1 先定等值面阈值，再谈好看
等值面是Liutex可视化最常用的方式。但阈值不能拍脑袋定，我的经验是先用Plot -> Contour把LiutexMag的分布范围看清楚，然后用最大值的一定比例来画等值面。实际操作中，我一般取LiutexMag最大值的10%~20%（例如最大值是150，阈值设在15到30之间），这样既能显示主涡核，又不会把次级弱涡全部淹没。如果是做论文图，建议给出2~3个不同阈值的对比图，审稿人最烦的就是“为什么选这个阈值”。
在Tecplot里生成等值面的路径是：Data -> Extract -> Iso-Surface，选LiutexMag，输入选定阈值。生成后可以切掉一半区域看内部结构，等值面显示用户通常会选Translucent，再用一个能区分正负旋转方向的颜色方案。
4.2 涡核线：直接沿Liutex向量追踪
Liutex向量天然给出了局部旋转轴方向r，这比从标量场提取涡核线要方便得多。在Tecplot里可以用Streamtrace功能，以Liutex向量分量(Lx, Ly, Lz)作为追踪向量，从等值面内部选若干个种子点，生成涡核线。这里有一个小技巧：种子点不要放在等值面的中心平面，而应该放在LiutexMag最大的截面上，用等值面和切平面交叉线来辅助定位种子点，得到的涡核线走向更稳定。
如果数据量比较大，建议先把LiutexMag场用插值或者仅提取局部区域的方式降采样再追踪，不然Tecplot的线积分速度会明显变慢。
4.3 高清图导出的关键设置：分辨率、矢量格式和字体
导出高清图片是投稿最刚需的操作。Tecplot 360导出图片时，File -> Export -> Image里面有几个选项必须确认：
- 输出格式：如果图形不复杂，优先选SVG或PDF等矢量格式，能保证放大后线条不糊；期刊要求PNG时再导出300 dpi的PNG。
- 位图分辨率：在PNG导出的
Resolution中选Custom，手动填300 dpi或更高。注意Tecplot里默认的Screen Resolution只有96 dpi，直接导出会被人吐槽“图片模糊”。 - 窗口尺寸：导图前把Tecplot窗口拖到你想要的宽高比，不要指望事后裁剪，Tecplot的导出跟当前视窗比例直接相关。
- 字体大小：导出前把坐标轴标签和色标字体统一调整到18磅以上，否则300 dpi下文字还会显得很小。
还有一个容易被忽略的点：Export面板里有个Supersample选项，可设置2×2或3×3超采样，能显著消除等值面边缘锯齿，但会大幅增加导出时间，适合最终定稿时用。
5. 我踩过的几个坑：边界梯度失真、复共轭混淆与多块数据错位
5.1 gradient在边界层和粗网格上的误差
Matlab的gradient函数在边界上用的是单侧差分，精度比内部中心差分低一阶。如果你分析的区域边界正好落在高剪切区，比如管流壁面或机翼表面，Liutex在壁面附近会出现虚假的高幅值。我试过忽略这个问题直接画图，结果等值面直接从壁面刮出一层伪涡壳，看起来像流动分离，其实是数值误差。
解决方法有两条：一是把计算域往无粘方向外扩几层虚拟网格，等梯度算完再截取物理区域；二是用更高精度的中心差分算子自己写梯度，比如用conv2或显式五点差分模板，避免边界单侧差分。
5.2 复共轭特征值接近实轴时的旋转方向误判
在剪切主导区域，速度梯度的复共轭特征对虚部很小，甚至因为数值误差出现虚部符号抖动。这时eig函数选出的“实特征值”可能不稳定，相邻网格点选出的特征向量会跳动。我通常会在符号一致性修正之前，先对λcr做一次滤波：把每个点与周围8个邻域的平均λcr做比较，超过3倍标准差就视为野点，用邻域中值替代。这个预处理对于涡量较弱区域的效果非常明显，等值面会干净很多。
5.3 ZONE方向写反导致的多块数据错位
多块数据最容易出现的问题不是算法，而是写文件时I/J/K顺序和Tecplot内部循环顺序不一致。如果某一块网格维度是[64, 128, 32]，但ZONE声明写成I=64, J=128, K=32，而输出循环还是按Matlab的行列页顺序，Tecplot会直接读错位置，画出完全不对的图。
我的排查技巧很笨但有效：导出后在文本文档里查看第一行ZONE设置，再抽查中间几行坐标数据，看X坐标是否单调递增。如果X方向坐标在文件中跳变，基本就是循环顺序的问题。
5.4 黑白印刷时的Liutex图识别度
最后分享一个和可视化审美相关的经验：期刊如果要求灰度图，Liutex等值面的正负旋转方向很容易在灰度印刷后失去区分度。我一般会切换Color Map为明暗对比较强的方案，比如从深灰到白色渐变，而不是默认蓝到红。另外可以叠加涡核线或矢量箭头，这样即使色标在灰度下看不清，结构特征也依然明确。
6. 我自己的使用心得：从“先画图后解释”到“先明确物理量再画图”
Liutex在Matlab+Tecplot里落地的整套流程，技术上并不算复杂，复杂的是你必须在每一步都知道自己在提取什么。我的体会是，不要在拿到CFD数据后立刻画Q准则或者Liutex等值面，而是先花半小时看速度梯度张量的统计特征，确定合理的阈值范围和旋转轴方向的一致性。数据的横截面探针、点云散点图这些“土办法”，往往能帮你省掉后期“图始终不干净”的数小时排错时间。做可视化不是炫技，是为了能说服包括审稿人在内的所有人：这个涡是真的存在，而且它的旋转方向、强度、位置都经得起检验。
本文还有配套的精品资源，点击获取
