简介：Liutex 方法作为第三代涡识别技术，常用于替代 Q 准则、λ2 等传统方法，这套资料面向流体力学研究人员与 CFD 工程师，提供从理论到代码实现的一站式参考。压缩包共 6 个文件、大小 6.67MB，其中 3 个 Matlab 脚本覆盖核心函数与数据写入流程，配套 PDF、PPTX 讲解阐明涡定义和实现思路，另有 eqn 步骤公式文件便于推导对照。已有 3582 人学习下载，适合希望快速掌握 Liutex 并衔接 Tecplot 可视化的中高级学习者。通过这份资料，可对照 Liutex 定义理解向量旋转与涡核提取要点，直接复用或修改 Matlab 脚本生成结果，再借助文档中的后处理步骤在 Tecplot 中显示涡结构，减少从公式到程序的搭建时间。资料包内同时提供 PDF 全文与 PPT 演示，便于从理论过渡到实践，尤其对需要做流动可视化报告或课题研究的读者很有帮助。

1. 为什么我卸掉了 Q 准则，改学 Liutex：一个新涡识别方法的前因后果

做完 CFD 后处理的人，大概率都被 Q 准则和 Lambda2 准则的阈值“玄学”折磨过——把阈值调得稍大，大涡结构稀碎；调得稍小，剪切层和边界层全被当真涡渲染出来，图一出来根本没法写进报告。Liutex 这个名字近两年在涡识别里越来越热，它最大的不同是直接刻画流场里真实的、相对刚体旋转的那一部分，而不是用某个阈值去等价成涡强度。配合 Matlab 算场、Tecplot 出图，这套实现流程完全可以复现到模拟项目 X 的后处理环节中。这篇文章就是我在给某跨平台系统做流场分析时，把 Liutex 从公式落实到可跑脚本的完整记录。适合手里有 CFD 数据、想摆脱阈值后处理烦恼的工程师。我们先从理论拆起，因为公式里的符号没吃透，后面代码里的参数就会变成纯黑匣子。

2. Liutex 矢量定义的物理拆解：从涡量到刚体旋转，核心到底是哪一步

2.1 为什么要单独定义“旋转”：涡量输运场里的水分太多

传统涡识别绕不开涡量 ω = ∇ × u。但涡量本身不区分剪切和刚体旋转，一个简单的平面库埃特流，涡量处处不为零，却没有一个真正的旋转核心。Q 准则把速度梯度张量拆成正对称和反对称部分，用应变率去抵消旋转贡献，本质上是个“能量平衡”的判断，所以它对阈值相当敏感。Omega 准则通过涡量平方和应变率平方之比把旋转占优的网格筛出来，鲁棒性好很多，但会丢掉局部旋转轴方向，这对湍流拟序结构分析很不友好。

Liutex 的定义思路完全不同：它认为流场中某点附近的运动可以唯一地分解成“刚体旋转部分”和“非旋转部分”，刚体旋转是真正的涡核特征。Liutex 是一个矢量 R，方向指向局部旋转轴，大小就是刚体旋转的角速度相关量，不再经过阈值筛选。换句话说，算完整个流场，每个网格点都会有一个明确的旋转矢量和模值，后续做等值面、做剖面统计都有直接的物理含义，不需要反复试阈值。这也是我后来把项目里的涡判据全部换成 Liutex 的根本原因。

2.2 局部旋转轴的确定：旋量代换法是怎么回事

Liutex 的计算难点在于求局部旋转轴方向 r。早期论文里给出了数学定义，但直接实现需要做特征值分解和约束优化。目前工程上最常用的是旋量代换法（curl replacement method），思路分两步：

第一步，先算涡量 ω，并归一化得到初始旋转轴方向 n = ω / |ω|。在大多数流动中，涡量方向和真实刚体旋转轴方向的偏差通常不大，可以作为迭代初值。

第二步，把速度梯度张量 G = ∇u 投影到以 n 为法向的平面和垂直平面上，提取“真正伴随旋转”的那部分。定义变换后的张量 G’，取其中包含旋转的二阶子张量，计算出局部角速度 λci（特征值虚部）。然后把旋量值更新为：

R = ω·r − sqrt( (ω·r)^2 − 4λci^2 )

如果 ω·r 不大于 2λci，说明这一点的变形以剪切为主，旋转贡献为零，直接赋 R = 0。这个公式的关键点在于：它把“涡量贡献”减去“剪切变形消耗”之后剩下的部分，才算作真实旋转。这么做的好处是数学上严格，物理上也不再依赖人工阈值。

2.3 参与计算的物理量：Rmag 与方向矢量 R 的分工

在实际写代码时，我会把 Liutex 结果拆成三个量来存储：Liutex 矢量本身 Rx、Ry、Rz，以及标量模值 Rmag。可视化时等值面尽量用 Rmag，因为它是标量，不会出现方向翻转导致显示断裂的问题；而分析涡轴走向时用 Rx/Ry/Rz，比如画涡管中心线，或者沿涡轴提取剖面。

有一点必须注意：Liutex 方向和涡量方向并不是完全一致的。涡量方向在剪切占主导的区域可能是个斜的指向，Liutex 方向因为做了二次修正，更贴近刚体旋转的真实轴。这也是为什么直接用涡量矢量做 iso-surface 时会看到很多“假核”，而换成 Liutex 矢量之后干净很多。后面 Matlab 代码里，我也会把这两个方向都导出来，方便大家自己对比差异。

3. Matlab 计算 Liutex：从三维速度场到 Rmag 的完整脚本与性能优化

3.1 数据准备：如何把 CFD 结果整理成 Matlab 能直接用矩阵

CFD 软件导出的数据有两种典型形态。非结构网格数据一般是一个点表和一个连接表，Matlab 处理起来要做散点插值；结构化网格数据是规则的三维数组，直接用索引差分就行。第一次尝试 Liutex 时，我建议先用结构化网格，把网格点坐标存成 X、Y、Z 三维数组，速度分量存成 U、V、W 三维数组，维度都是 Nx × Ny × Nz。如果你的原始数据是非结构的，先用 griddata 或 scatteredInterpolant 插值到规则网格上，这一步虽然损失一点精度，但能让你把算法调通。

常见的坑是网格间距不一致。很多程序只会输出 x/y/z 坐标，不会直接给出 dx、dy、dz。我一般会检查坐标数组，如果坐标是线性分布就直接算间距标量，如果是非均匀网格，就生成三个与网格同尺寸的间距数组。中心差分时，非均匀间距的差分系数不能简单套均匀网格公式，否则梯度场会有锯齿。下面的代码先按均匀网格处理，非均匀情况我在后面避坑章节里单独说明。

3.2 速度梯度与涡量计算：用中心差分避开内置方法的限制

Matlab 自带 curl 函数，但它是针对二维平面场设计的，三维情况下不能省事直接调。同时它也把方向定义为正切于坐标平面的旋转，不符合流体力学里的右手定则标注。所以我习惯自己写梯度函数，用中心差分：

function liu = compute_liutex(U, V, W, dx, dy, dz) % 输入：U/V/W 三维速度分量数组 % dx/dy/dz 网格间距标量或与 U 同尺寸的数组 % 输出：Rx/Ry/Rz 分量与 Rmag 绝对值 [Nx, Ny, Nz] = size(U); % 预分配空间，避免循环里动态扩维 Rx = zeros(Nx, Ny, Nz); Ry = zeros(Nx, Ny, Nz); Rz = zeros(Nx, Ny, Nz); Rmag = zeros(Nx, Ny, Nz); % 边界用前向/后向差分，内部用中心差分 for k = 2:Nz-1 for j = 2:Ny-1 for i = 2:Nx-1 duidx = (U(i+1,j,k) - U(i-1,j,k)) / (2*dx); duidy = (U(i,j+1,k) - U(i,j-1,k)) / (2*dy); duidz = (U(i,j,k+1) - U(i,j,k-1)) / (2*dz); dvidx = (V(i+1,j,k) - V(i-1,j,k)) / (2*dx); dvidy = (V(i,j+1,k) - V(i,j-1,k)) / (2*dy); dvicz = (V(i,j,k+1) - V(i,j,k-1)) / (2*dz); dwidx = (W(i+1,j,k) - W(i-1,j,k)) / (2*dx); dwidy = (W(i,j+1,k) - W(i,j-1,k)) / (2*dy); dwidz = (W(i,j,k+1) - W(i,j,k-1)) / (2*dz); % 构造速度梯度张量 G G = [duidx, duidy, duidz; dvidx, dvidy, dvicz; dwidx, dwidy, dwidz]; % 涡量 omega = curl(V) wx = (V(i,j,k+1)-V(i,j,k-1))/(2*dz) - (W(i,j+1,k)-W(i,j-1,k))/(2*dy); wy = (W(i+1,j,k)-W(i-1,j,k))/(2*dx) - (U(i,j,k+1)-U(i,j-1,k))/(2*dz); wz = (U(i,j+1,k)-U(i,j-1,k))/(2*dy) - (V(i+1,j,k)-V(i-1,j,k))/(2*dx); omega_norm = sqrt(wx^2 + wy^2 + wz^2); % 如果涡量为零，则这一片没有旋转 if omega_norm < 1e-10 continue; end % 初始旋转轴方向取涡量归一化方向 nvec = [wx, wy, wz] / omega_norm; [R, ~, ~] = liutex_from_gradient(G, nvec); Rx(i,j,k) = R(1); Ry(i,j,k) = R(2); Rz(i,j,k) = R(3); end end end Rmag = sqrt(Rx.^2 + Ry.^2 + Rz.^2); end

这块代码逻辑很直白：对每个内部网格点做中心差分求速度梯度张量，再求涡量归一化方向 nvec，最后把梯度张量和方向交给核心子函数liutex_from_gradient去算。边界层暂时跳过，否则索引会越界。

需要注意：如果 dx、dy、dz 是标量，代码没问题；如果是非均匀网格的数组，这里要改成(U(i+1,j,k)-U(i-1,j,k)) / (x(i+1)-x(i-1))这种索引式写法。另外，涡量为零的点直接 continue 掉，省去后面无意义的矩阵运算。

3.3 旋量代换法核心子函数：解析计算 λci，避免硬扛特征值分解

子函数里最重要的一步是把速度梯度张量沿 nvec 方向做坐标变换。我的做法是构造一个旋转矩阵 Q，让 nvec 对应新坐标系的 x 轴，然后把 G 相似变换到新坐标系下。变换后 G’ 的 2x2 子块包含了垂直于旋转轴平面上的变形信息。此时可以直接用解析公式求 λci：

λci = sqrt( max(0, 4*det - trace^2) ) / 2

其中 det 和 trace 都是那个 2x2 子块的行列式与迹。解析公式比逐点调用 eig 快两个数量级，大阵列仿真的性能瓶颈一般就在这一步。

function [Rvec, Rscalar, flag] = liutex_from_gradient(G, nvec) % 输入：速度梯度张量 G (3x3)，单位旋转轴方向 nvec % 输出：Liutex 矢量 Rvec，标量 Rscalar，以及是否有效标志 flag % 构造正交基，使 nvec 为新 x 轴 e1 = [1 0 0]; if abs(dot(nvec, e1)) > 0.9 e1 = [0 1 0]; % 如果 nvec 与原 x 轴接近，改用 y 轴为辅助 end e2 = cross(nvec, e1); e2 = e2 / norm(e2); e1 = cross(e2, nvec); % 重新正交化 % Q 的行是新基向量，Q*G*Q' 把梯度张量投影到新坐标系 Q = [nvec; e1; e2]; G_rot = Q * G * Q'; % 提取垂直于旋转轴的 2x2 子块 A = G_rot(2:3, 2:3); trA = A(1,1) + A(2,2); detA = A(1,1)*A(2,2) - A(1,2)*A(2,1); lambda_ci = sqrt(max(0, 4*detA - trA^2)) / 2; % 涡量在轴向的投影 omega_dot_r = G(3,2) - G(2,3); % 涡量沿 nvec 方向的分量 % 旋量代换法核心公式 inner = omega_dot_r^2 - 4*lambda_ci^2; if inner < 0 Rscalar = 0; flag = 0; % 剪切主导，不认为存在刚体旋转 else Rscalar = omega_dot_r - sqrt(inner); if Rscalar < 0 Rscalar = 0; flag = 0; else flag = 1; end end Rvec = Rscalar * nvec; end

从代码可以看到，旋转坐标系的构造顺序很关键：先用 nvec 叉乘辅助向量得到 e2，再用 e2 叉乘 nvec 得到 e1，这样可以保证三个基向量构成右手系。如果顺序颠倒，旋转轴方向会整体翻转，导出的 Liutex 矢量全都会指向反方向。我在最初实现时就在这里翻过车，图里涡轴方向在相邻网格间来回跳，看起来像信号紊乱，其实只是坐标系构造顺序问题。

lambda_ci这里用的是二维子块的解析公式，完全绕开了eig函数的复杂分支。4*detA - trA^2如果不加 max(0,...) 保护，在剪切占主导的区域会出现负值，sqrt 直接给 NaN，后面无论怎么可视化都是白屏，这是最常见的内存级“暗坑”。

3.4 性能瓶颈与向量化方案：别在三维循环里逐点计算

上面这个版本适合先跑通逻辑。如果网格规模是 500×300×200，逐点三重循环加坐标变换，Matlab 里可能要跑十几分钟。常见做法是把循环体改成逐列或逐面操作。我发现直接去掉中间层循环，把 j 方向向量化，能稳定提速 5~8 倍：

% 以 k 维度最外层循环，j 和 i 采用数组切片 for k = 2:Nz-1 dudx = (U(3:Nx,j,k) - U(1:Nx-2,j,k)) / (2*dx); dudy = (U(2:Nx-1,3:Ny,k) - U(2:Nx-1,1:Ny-2,k)) / (2*dy); % 以二维矩阵块方式计算 G ... end

这个版本的写起来稍显繁琐，但在数据量大的时候非常值得。我一般先把三维数组降成若干二维面，逐面调用一个矩阵化函数，避免 Matlab 的 JIT 编译器在深层循环里矫情。如果要处理大规模 DNS 或者 LES 数据，甚至可以把梯度计算部分改成 mex 函数，或者用 GPU 数组。不过第一次上手，先保证和理论结果对上，再去折腾性能。

4. 从 Matlab 到 Tecplot：把 Liutex 结果导出并渲染出能放进论文的图

4.1 数据导出格式选择：ASCII dat 是最稳妥的互通方案

Tecplot 读取数据的方式很多，但最不挑版本的就是 ASCII 的 .dat 文件。它本质上是文本格式，Matlab 用 fprintf 写出来，Tecplot 直接 File -> Load Data Files 就能读。如果数据量特别大，比如超过 2000 万网格点，文本文件确实会非常占磁盘，可以考虑 Tecplot 的二进制 plt 格式，但初次实验完全没有必要。

为了后续切片和等值面方便，我的导出脚本会把 X、Y、Z、U、V、W、Rx、Ry、Rz、Rmag 这十个变量都写进同一个文件。这样 Tecplot 里一个 zone 就包含全部数据，切换显示变量比反复导多次文件省事很多。坐标单位要和速度单位保持一致，这个很多人会忽略，导致后期算 Q 准则的公式时量纲对不上。

function export_to_tecplot(filename, X, Y, Z, U, V, W, Rx, Ry, Rz, Rmag) % 导出结构化网格 Liutex 结果到 Tecplot ASCII 文件 [Nx, Ny, Nz] = size(X); fileID = fopen(filename, 'w'); fprintf(fileID, 'TITLE="Liutex simulation case"\n'); fprintf(fileID, 'VARIABLES="X","Y","Z","U","V","W","RX","RY","RZ","RMAG"\n'); fprintf(fileID, 'ZONE I=%d, J=%d, K=%d, ZONETYPE=ORDERED, DATAPACKING=POINT\n', Nx, Ny, Nz); for k = 1:Nz for j = 1:Ny for i = 1:Nx fprintf(fileID, '%e %e %e %e %e %e %e %e %e %e\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), ... Rx(i,j,k), Ry(i,j,k), Rz(i,j,k), Rmag(i,j,k)); end end end fclose(fileID); end

这段写入逻辑是逐个网格点输出一行。注意 Tecplot 的ZONETYPE=ORDERED必须和 Matlab 数组的索引顺序对应，不然一进入 Tecplot 就会看到网格拓扑完全乱掉。DATAPACKING=POINT 表示一行一个数据点，清晰且不容易错；如果选中 BLOCK 模式，需要把同一变量所有数据一次性输出，文本行会非常长，不利于人工检查第一行数据是否合理。

导完以后我还会用文本编辑器打开文件头几行，确认第一行坐标不是负数、速度没有出现 NaN，再做后续操作。这一步能提前挡住大量渲染时的“莫名其妙白屏”问题。

4.2 Tecplot 可视化的三个关键步骤：切片、等值面、轴向着色

拿到数据后，我会在 Tecplot 里按固定套路操作。

第一步，做剖面。进入 Plot -> Slice，选取 X-Y 平面切一刀，把 Z 方向速度为 0 的剖面显示出来。变量选择 Rmag，Zones 限制在当前唯一的一个 zone 上。切片的好处是能快速确认涡核位置分布是否符合物理直觉，比如绕圆柱体的卡门涡街，切片上应该能看到周期性的交替脱落。

第二步，做等值面。Tecplot 的等值面默认基于某变量，这里选 Rmag。和传统 Q 准则等值面不同，Rmag 的数值有明确的“刚体旋转强度”含义，所以我们可以用同一套阈值处理不同工况的数据，不再需要每张图重新调。我一般先取全局 Rmag 最大值的 20% 作为初始值，然后上下微调 5%，找到涡核清晰且连接体完整的水平。

第三步，给等值面上色。这个上色变量可选的很多，常见做法是选 velocity magnitude，或者根据流场方向选某个速度分量。如果关心流向变化，我习惯用 U 分量着色；如果关心旋转强度，干脆用 Rmag 自己给自己着色，得到的是从红到蓝的“涡强度云图”。Tecplot 的色标范围默认自动适应，但为了论文统一，必须手动设置最大最小值。

4.3 用 Tecplot 宏批量出图：一键生成多个截面的动画脚本

如果想看三维涡结构在不同截面上的演化，手动操作几十个切片会累死人。Tecplot 支持宏脚本，可以把操作录下来重放。我常用的方式是用 Tecplot 的宏录制功能先手动操作一次，再在文本编辑器里修改截面位置参数，得到类似下面的模板：

#!MC 1410 $!REDRAWALL $!CREATEEQUATION SOLVEEQUATION = 'XSLICE = 0.1' $!SETACTIVESTRAND = 1 $!EXTENDSLICES

注意#!MC 1410是 Tecplot 宏文件的固定头部格式，版本不同的话后缀会有区别。我不推荐直接手写宏，建议先在界面上录一遍，再在生成的脚本里改截面坐标和文件输出路径。录制方式是 Macros -> Record Macro，做完切片后选保存，Tecplot 会把整个交互过程转化为可编辑的脚本。批量出图时，把循环变量改成数组，就能直接跑出一系列 PNG 帧，后期用视频工具拼成动画。

Tecplot 切片面板里的“Extend To Domain”选项也有讲究。如果切面没有延伸到几何域的所有方向，最终图像外圈会留下一圈空白；我一般保持默认不勾选，让切面和整个计算域重叠，这样不同工况的图在视觉上是同等范围，对比起来不会因为边界差异产生误导。

5. 实战避坑：Liutex 计算与可视化中让我往返修改的五类问题

5.1 中心差分边界处理导致涡核“断裂”

现象：等值面在大涡核中间出现规则性的缺口，但涡量切片显示流场明明是连续旋转的。缺口位置恰好位于i=2或j=Ny-1这类边界附近，越靠近外侧越明显。

原因：核心算法里我把边界点跳过了。当计算域边缘的导数值只由一侧决定时，我直接用前向差分，而内部点用中心差分，两种差分格式的截断误差不同，造成 Liutex 幅值在边界附近突然掉下来，等值面迎合一截就断。

解决：边界层附近建议至少保留三层网格做差分缓冲。如果计算域本身的涡核靠近边界，就用多项式拟合求导，而不是简单的前向/后向差分。我后来在代码里写了边界处理分支：当i==1时用(-3U(1,j,k)+4U(2,j,k)-U(3,j,k))/(2dx)，这个二阶精度单侧差分公式比普通前向差分平滑很多，边界处的 Rmag 不会出现断崖。

5.2 根号内出现负数导致大量 NaN，整个云图变白屏

现象：Matlab 脚本跑完，检查结果矩阵里全是 NaN，Tecplot 导入后画面一片空白，连等值面都没有。检查代码发现sqrt(omega_dot_r^2 - 4*lambda_ci^2)经常报“复数的 sqrt 不推荐”警告。

原因：当局部变形张量的剪切贡献超过了涡量的有效旋转贡献时，根号内确实为负。这在理论上正好说明该点不是合适的旋转核心。不加保护的 sqrt 会把实数域扩展到复数域，得到的虚数无法参与后续绘图。

解决：在核心函数里对根号内做max(0, ...)钳位。如果结果为负，直接置 Rscalar=0，这既符合物理判断，也让代码在数值上变得健壮。我那版代码里inner < 0分支已处理了这层逻辑，后续所有网格点都没有再出现 NaN。

5.3 非结构化网格直接套用三维数组索引，结果完全错乱

现象：从某软件导出的 CFD 数据是三角形/四面体网格，点数和网格数不一致。直接把散点坐标 reshape 成三维数组，Tecplot 里看网格拓扑完全是乱线，等值面位置也完全对不上。

原因：结构化网格的索引连续性在非结构网格里不存在。我虽然用 griddata 做了插值，但插值时的散点坐标顺序没有和原始网格的局部编号对齐，导致插值结果其实是“张冠李戴”。

解决：插值前先检查原始点云的空间包围盒。把坐标中心平移到原点，再统一单位。我习惯先用scatteredInterpolant建立速度场的插值器，再在规则网格的坐标上做评估，最后导出之前再画一张快速切片对比原始数据，确认没有局部错位。非结构网格过渡到规则网格后，Liutex 计算精度取决于插值方法，线性插值一般够用，但要注意网格分辨率不能比原始网格粗太多。

5.4 Liutex 方向矢量在涡心附近“打架”，流线图出现反向跳跃

现象：渲染涡轴流线时，相邻两条流线在涡心区域突然反向，颜色也由红变蓝。单独看每个点 Rmag 幅值似乎正常，就是方向不稳定。

原因：Liutex 矢量方向在局部旋转轴穿过零涡量线时，会因为数值误差出现符号翻转。虽然理论上方向连续，但数值计算中涡量方向接近零时，归一化方向对舍入误差极其敏感，稍微一点不对称扰动就能让方向跳到反方向。

解决：做流线显示时不要直接画矢量头，而是绘制基于Rmag 等值面 + Rx/Ry/Rz 分量投影的组合。比如 Rmag 等值面确定涡核位置，再用速度矢量在场内做流线，加上的箭头的疏密靠 Rmag 调制，就不会出现方向跳变导致的视觉撕裂。如果必须直接显示 Liutex 矢量场，可以先用高斯滤波把方向场平滑一遍。

5.5 大网格数据内存爆掉，运行到一半 Matlab 直接崩溃

现象：三维数组分别是 1024×512×512 时，代码跑到三分之一就提示内存不足，系统物理内存使用了 90% 以上。

原因：预分配 Rx、Ry、Rz 和 Rmag 四个双精度矩阵，每个都是 1024×512×512×8 字节约 2GB，四个就是 8GB。加原始速度场和梯度场，整机 32GB 内存也扛不住。

解决：Rmag 和矢量的存储可以降精度为 single 类型，精度损失在可视化尺度上几乎不可见。更彻底的做法是分块计算：按 k 方向切出若干子块，每个子块单独调 compute_liutex 函数，算完直接写入磁盘或 Tecio 导出文件，不保留全局结果矩阵。我通常把 k 循环切成 64 层一块，循环处理完所有切片后再整合。实测内存占用能降到 5GB 以内，运行时长只多了约 8% 的边界重复计算。

6. 验证与更进阶的玩法：用 Taylor-Green 涡检验代码，并把 Liutex 方向信息用于涡轴追踪

拿到一套新算法代码，第一件事不是去算大工程，而是先做数学验证。我会构造一个已知解析解的 Taylor-Green 涡场：速度分量满足 u = sin(x)cos(y)cos(z)，v = -cos(x)sin(y)cos(z)，w = 0。理论上 Liutex 幅值在涡心位置应该等于对应角速度的两倍，并且方向严格对齐 z 轴。把这个理论场喂给上面的 Matlab 代码，对比 Rmag 解析解和数值解，最大误差控制在 1% 以内，就说明代码整体是可靠的。

这个验证步骤看起来简单但非常重要。我在给模拟项目 X 部署到新机器后，重新运行时发现数值结果比之前差了 3%，排查半天发现是编译器的浮点优化选项变了，导致差分精度下降。从那以后，我每次修改任何参数或切换环境，都会先跑一遍 Taylor-Green 验证案例，给自己一个“后悔药”级别的安全网，再开展大规模分析。

验证通过后，还可以进一步用 Liutex 矢量方向追踪涡轴。做法是：从 Rmag 等值面上取最高幅值点作为种子，然后沿局部 Liutex 方向做步进，每走一步在该点附近重新计算 Rmag 和方向，持续积分得到涡轴中心线。这个方法比传统压力极小值法稳定，尤其在高剪切环境下，压力极值点常常不在真实涡心。如果是研究边界层转捩，还能用 Liutex 矢量方向做条件平均，观察不同涡片段的旋转轴统计分布，这是 Q 准则完全给不了的信息。

数据显示管路的最终习惯是：我用 Matlab 脚本把 Nd 帧动画全部算好，Tecplot 里用一个循环宏批量导出 PNG，再用工具把 PNG 合成为动画。动画里能看到涡核从产生、拉伸到破裂全过程，而且因为 Liutex 已经在物理上剔除了剪切，不会有剪切层导致的视觉噪音。希望这份记录能帮到准备在 Tecplot 里把涡识别做得更干净的人。

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