简介这份资源围绕Liutex涡识别方法展开面向流体力学、湍流研究及CFD后处理方向的学习者与科研人员帮助理解第三代涡识别方法的数学定义并借助Matlab与Tecplot完成从计算到可视化的完整流程。压缩包共6个文件约6.67MB包含3个m脚本、1个pptx演示文稿、1个eqn公式文件与1个pdf文献脚本负责Liutex相关量的计算与写出公式文件便于在文档中复用推导演示与文献则补充理论背景。资源中提供了Liutex函数实现、旋转量提取与结果输出等核心代码读者可据此复现涡结构识别过程并将其接入自有流场数据配合Tecplot完成涡结构可视化与对比分析。目前已有3576人学习下载适合希望快速上手Liutex方法、减少重复推导与编码工作的研究者参考使用。1. 从涡旋识别到可视化落地Liutex 方法在 Matlab 与 Tecplot 中的实现路径做湍流或者旋涡分离的朋友大概率都经历过这种场景算完一个叶轮机械或者钝体绕流的瞬态场想从速度场里把涡结构拎出来Q 准则、λ2 准则、Δ 准则轮番上阵结果要么阈值调得人崩溃要么涡核边界糊成一团。Liutex 方法这两年被提得越来越多核心原因就是它给出了一个相对严格的刚体旋转分量定义不再依赖经验阈值去“猜”涡在哪。这套资源围绕的就是怎么在 Matlab 里把 Liutex 算出来再顺到 Tecplot 里做场可视化。它适合已经能跑 CFD、手里有速度场数据、但卡在涡识别后处理这一步的从业者也适合想从 Q 准则迁移到 Liutex 但不知道从哪下手的同学。Matlab 负责矩阵运算和批量处理Tecplot 负责最终出图和动画中间的数据接口和坐标系对齐是整条链路里最容易翻车的地方。2. Liutex 的数学内核从速度梯度张量到刚体旋转轴2.1 为什么 Q 准则不够用局部旋转与剪切的对立Q 准则的定义是速度梯度张量的第二不变量为正物理上对应旋转率大于应变率。问题在于剪切层里 Q 值也可以很大导致你把剪切误判成涡。λ2 准则稍微好一点看的是压力张量 Hessian 的次特征值但它本质上还是在找“压力极小面”对可压缩或者强旋流的适应性有限。Liutex 的思路不一样它试图把局部速度梯度分解成一个刚体旋转部分和一个非刚体部分然后直接给出旋转轴的向量和旋转强度。这个分解不是唯一的但 Liu 等人给出的定义在数学上满足伽利略不变性而且旋转强度是量纲一致的。实际用下来最大的好处是你不用再对着阈值反复试Liutex 的等值面在涡核附近梯度更陡边界更干净。2.2 速度梯度张量的实特征值与复特征值判定Liutex 的计算起点是速度梯度张量 ∇v一个 3×3 矩阵。先求它的特征值如果三个特征值全是实数说明局部没有旋转Liutex 为零如果有一个实特征值和一对共轭复特征值实特征值对应的特征向量就是旋转轴方向。这一步在 Matlab 里用 eig 就能做但要注意数据排列顺序。很多 CFD 导出的是扁平化的速度场需要先 reshape 成 [nx, ny, nz, 3] 再逐点算梯度。梯度计算建议用中心差分边界用单侧差分不要直接用 gradient 函数因为非均匀网格下 gradient 默认假设间距为 1会引入系统性偏差。% 假设 u, v, w 是 nx*ny*nz 的速度分量x, y, z 是对应坐标 [nx, ny, nz] size(u); % 预分配速度梯度张量 grad_v zeros(nx, ny, nz, 3, 3); % 中心差分计算梯度循环维度可改为向量化 for i 2:nx-1 for j 2:ny-1 for k 2:nz-1 dx x(i1,j,k) - x(i-1,j,k); dy y(i,j1,k) - y(i,j-1,k); dz z(i,j,k1) - z(i,j,k-1); grad_v(i,j,k,1,1) (u(i1,j,k) - u(i-1,j,k)) / dx; grad_v(i,j,k,1,2) (u(i,j1,k) - u(i,j-1,k)) / dy; grad_v(i,j,k,1,3) (u(i,j,k1) - u(i,j,k-1)) / dz; grad_v(i,j,k,2,1) (v(i1,j,k) - v(i-1,j,k)) / dx; grad_v(i,j,k,2,2) (v(i,j1,k) - v(i,j-1,k)) / dy; grad_v(i,j,k,2,3) (v(i,j,k1) - v(i,j,k-1)) / dz; grad_v(i,j,k,3,1) (w(i1,j,k) - w(i-1,j,k)) / dx; grad_v(i,j,k,3,2) (w(i,j1,k) - w(i,j-1,k)) / dy; grad_v(i,j,k,3,3) (w(i,j,k1) - w(i,j,k-1)) / dz; end end end这段代码的逻辑很直接对每个内部网格点用相邻两点的速度差除以坐标差得到偏导数。参数说明上dx、dy、dz 是局部网格间距非均匀网格必须逐点算不能用一个全局平均值代替。边界点单独处理通常用前向或后向差分。如果网格量级在百万以上三重循环会非常慢建议用 permute 和 diff 做向量化或者干脆写 mex。我一般会先把梯度张量存成稀疏格式因为大部分点的 Liutex 值接近零稠密存储浪费内存。2.3 复特征值对应的旋转轴提取拿到 ∇v 之后逐点求特征值。Matlab 的 eig 对 3×3 矩阵返回一个 3×1 复向量。判断是否有复特征值如果任意两个特征值的虚部绝对值大于某个小量比如 1e-12就认为存在旋转。实特征值对应的特征向量就是旋转轴。注意eig 返回的特征向量可能符号相反需要统一方向否则 Tecplot 里画出来的涡轴会跳变。常见做法是让旋转轴与当地涡量方向点积为正如果为负就取反。liutex_mag zeros(nx, ny, nz); liutex_axis zeros(nx, ny, nz, 3); for i 1:nx for j 1:ny for k 1:nz A squeeze(grad_v(i,j,k,:,:)); [V, D] eig(A); eig_vals diag(D); % 判断是否存在复特征值 if any(abs(imag(eig_vals)) 1e-12) % 找实特征值对应的特征向量 [~, idx] min(abs(imag(eig_vals))); axis_vec real(V(:, idx)); axis_vec axis_vec / norm(axis_vec); % 与涡量对齐 vort [grad_v(i,j,k,3,2)-grad_v(i,j,k,2,3), ... grad_v(i,j,k,1,3)-grad_v(i,j,k,3,1), ... grad_v(i,j,k,2,1)-grad_v(i,j,k,1,2)]; if dot(axis_vec, vort) 0 axis_vec -axis_vec; end liutex_axis(i,j,k,:) axis_vec; % Liutex 强度计算这里用简化形式 liutex_mag(i,j,k) abs(imag(eig_vals(idx))); end end end end这里 liutex_mag 只是取了虚部的绝对值作为强度的一个代理严格定义下 Liutex 强度是旋转轴方向上的速度梯度分量减去其转置分量。实际写的时候建议按原始文献的公式来我这里为了展示流程做了简化。参数上1e-12 这个容差要根据数据精度调双精度一般够用单精度数据要放大到 1e-6 左右。liutex_axis 存成四维数组后面导出 Tecplot 的时候直接对应三个分量。3. Matlab 到 Tecplot 的数据接口格式选择与坐标系对齐3.1 Tecplot 能吃的几种格式PLT、DAT 与 CGNSTecplot 支持的原生格式里二进制 PLT 读写最快但 Matlab 没有内置写 PLT 的函数得靠第三方 mex 或者自己拼二进制头。DAT 是 ASCII 格式写起来简单但文件体积大百万网格点轻松上 G。CGNS 是 CFD 通用格式Tecplot 和大部分 CFD 软件都认但 Matlab 写 CGNS 需要额外工具箱。我一般用 DAT 做小规模调试正式出图用 PLT。如果你们组有 Tecplot 的 Chorus 或者 360 版本也可以走 Python 的 pytecplot 接口但那就绕开 Matlab 了。这里只讲纯 Matlab 写 DAT 的路子因为最可控。3.2 用 Matlab 拼 Tecplot ASCII 数据文件Tecplot 的 DAT 文件结构分三块标题区、变量列表、数据区。变量列表里要声明所有你要画的量比如 X、Y、Z、U、V、W、Liutex、Axis_X、Axis_Y、Axis_Z。数据区按变量顺序逐点排列。对于结构网格可以用 I、J、K 声明维度然后按 i 最快、j 次之、k 最慢的顺序写。下面是一个写 DAT 的示例。fid fopen(liutex_output.dat, w); fprintf(fid, TITLE Liutex Field\n); fprintf(fid, VARIABLES X Y Z U V W Liutex Ax Ay Az\n); fprintf(fid, ZONE TFlowField, I%d, J%d, K%d, DATAPACKINGPOINT\n, nx, ny, nz); for k 1:nz for j 1:ny for i 1:nx fprintf(fid, %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), ... liutex_mag(i,j,k), ... liutex_axis(i,j,k,1), liutex_axis(i,j,k,2), liutex_axis(i,j,k,3)); end end end fclose(fid);逻辑说明先写标题和变量名Tecplot 对变量名大小写不敏感但建议统一。ZONE 行里的 I、J、K 告诉 Tecplot 这是结构网格DATAPACKINGPOINT 表示按点写。循环顺序必须是 k 外层、j 中层、i 内层因为 Tecplot 默认 i 变化最快。参数上%e 是科学计数法精度够用如果数据范围跨好几个量级可以改用 %g。文件写完后用 Tecplot 打开如果报“数据点数量不匹配”九成是循环顺序反了。3.3 非结构网格的导出策略如果你的网格是非结构的比如四面体或者多面体DAT 文件里就不能用 I、J、K 了得改成 ZONE 里声明 NODES 和 ELEMENTS然后先写节点坐标和变量值再写单元连接关系。Matlab 里通常有网格的节点坐标和单元索引矩阵直接按格式拼就行。注意 Tecplot 的单元索引是从 1 开始的Matlab 也是从 1 开始这点不用转换。但如果你中间用过其他工具转过格式要确认索引基准没变。非结构网格的 Liutex 轴可视化建议用切片或者流线不要直接画等值面因为单元形状不规则等值面容易破。4. Tecplot 里的 Liutex 可视化等值面、涡轴与动画4.1 用 Liutex 幅值做等值面阈值怎么选数据导入 Tecplot 后第一件事是画 Liutex 的等值面。在 Contour 里选 Liutex 变量然后切到 Isosurface 模式。阈值选取有个经验先看 Liutex 的体渲染或者切片云图找到涡核区域的典型值然后取该值的 30% 到 50% 作为等值面阈值。Liutex 的好处是涡核内部值高、外部衰减快所以阈值窗口比 Q 准则窄不用反复试。如果等值面出现碎片说明阈值偏低把剪切区的噪声也包进来了如果涡结构断成几截说明阈值偏高把涡核弱连接处切掉了。我一般会同时开两个等值面一个高阈值看涡核一个低阈值看涡的包裹范围。4.2 涡轴向量场的显示用 Vector 还是 StreamtraceLiutex 轴是单位向量直接画 Vector 图会密密麻麻而且箭头长度没有物理意义。更好的做法是用 Streamtrace 沿着涡轴方向积分或者用 Tecplot 的 Vortex Core 功能线。具体操作在 Plot 里选 Vector变量选 Ax、Ay、Az然后调整箭头间距和长度。如果涡轴方向跳变严重检查 Matlab 里有没有做方向对齐。另一个技巧是用 Volume 渲染把 Liutex 幅值映射成透明度涡轴用 Line 画这样既有强度信息又有方向信息。Tecplot 的 Volume 渲染对显卡有要求如果卡顿就改用切片。4.3 瞬态动画的帧序列导出瞬态数据通常是一个时间序列一个文件Tecplot 可以用 Time 菜单里的 Animate 做帧切换。前提是每个 DAT 文件里的变量名和网格维度完全一致否则 Tecplot 会认为是不同数据集。导出动画时选 AVI 或者 GIF帧率一般 10 到 15 帧每秒就够了。如果文件太多可以用 Tecplot 的宏命令批量加载和导出。宏文件里用 $!ReadDataSet 读文件$!Export 导出图片循环变量用 $!Loop。注意宏里的路径要用双反斜杠或者正斜杠Windows 下单反斜杠会被当转义符。5. 避坑与排查Liutex 计算和可视化里的五个血泪教训5.1 现象Liutex 全场为零但流场明明有涡原因速度梯度张量的特征值全是实数说明局部没有复特征值。常见于涡核被数值耗散抹平或者网格分辨率不够涡的旋转特征在离散尺度上退化了。解决先检查网格在涡核附近的尺寸如果当地网格尺度大于涡核半径Liutex 必然算不出来。加密网格或者用高阶插值重构速度场。另一个可能是速度数据本身有问题比如从 CFD 导出时做了太多平滑。5.2 现象Tecplot 打开 DAT 文件报“变量数量不匹配”原因VARIABLES 行声明的变量个数和实际每行写的数据个数不一致。常见于变量列表里多写了一个空格或者少写了一个引号。解决用文本编辑器打开 DAT 文件数一下 VARIABLES 行里引号对的数量再数数据行里数值的个数。Tecplot 对格式很敏感多一个逗号都可能报错。建议在 Matlab 里用 fprintf 拼字符串时把变量名和格式串分开维护改一个地方就同步改另一个。5.3 现象涡轴方向在相邻点之间来回翻转原因eig 返回的特征向量符号是任意的Matlab 不保证连续性。解决在计算每个点的轴向量后与上一个点的轴向量做点积如果为负就取反。更稳健的做法是与当地涡量方向对齐因为涡量方向在涡核附近是连续的。如果涡量本身也跳变那就用全局参考方向比如来流方向但只适用于涡轴大致一致的场景。5.4 现象Tecplot 里等值面出现大量碎片像噪声原因Liutex 在剪切层里虽然比 Q 准则干净但仍有小值残留。如果阈值取得太低这些残留就被包进来了。解决先对 Liutex 场做一次空间平滑用 3×3×3 的均值滤波但注意平滑会削弱涡核峰值。另一个办法是在 Matlab 里设一个下限低于某个值的 Liutex 直接置零这个下限取全场最大值的 1% 到 5%。置零后等值面会干净很多。5.5 现象瞬态动画里涡结构闪烁不连续原因每帧的 Liutex 阈值是独立选的或者每帧的坐标轴范围变了。解决在 Tecplot 里锁定等值面阈值和视图范围所有帧用同一套参数。如果数据本身有时间步之间的数值振荡可以在 Matlab 里对时间序列做低通滤波。另外瞬态数据的 DAT 文件命名要有规律比如 liutex_0001.dat、liutex_0002.datTecplot 的 Time 菜单才能自动识别序列。6. 进阶技巧用 Matlab 批量算 Liutex 并直接驱动 Tecplot 宏当你需要处理几十个工况或者几百个时间步的时候手动在 Tecplot 里点来点去就不现实了。我一般会在 Matlab 里写一个主脚本循环读速度场、算 Liutex、写 DAT然后调用 Tecplot 的宏做批量出图。Tecplot 宏可以用命令行调用Windows 下是 tecplot.exe -b -p macro.mcrLinux 下类似。宏文件里用 $!ReadDataSet 读 DAT$!RedrawAll 刷新$!Export 导出图片。下面是一个宏模板。# Tecplot 宏模板批量导出 Liutex 等值面图片 $!ReadDataSet liutex_output.dat $!CreateContourGroup [1] $!ContourGroup [1] ContourType Isosurface $!ContourGroup [1] ContourVar 7 $!ContourGroup [1] IsosurfaceValue 0.35 $!RedrawAll $!ExportSetup ExportFormat PNG $!ExportSetup ImageWidth 1920 $!ExportSetup ImageHeight 1080 $!Export liutex_iso.png这个宏的逻辑是读数据建等值面组指定变量索引为 7对应 Liutex设等值面值为 0.35刷新后导出 PNG。参数上ContourVar 的索引要和 VARIABLES 行里的顺序一致从 1 开始数。IsosurfaceValue 根据你的数据范围调可以先跑一个工况看云图再定。ExportSetup 里的分辨率按期刊要求设一般 300 dpi 对应 1920 宽就够。如果要做动画把 ReadDataSet 和 Export 放进 $!Loop 里循环变量用文件名后缀。还有一个技巧是直接在 Matlab 里生成 Tecplot 的布局文件.lay这样连宏都不用写。布局文件是 XML 格式里面定义了帧、等值面、切片、流线等所有可视化元素。Matlab 可以用 xmlwrite 或者直接 fprintf 拼。好处是一次生成以后每次只换数据文件就行。缺点是布局文件的 schema 在不同 Tecplot 版本之间有差异换版本可能要微调。我一般会在目标版本上手动存一个模板然后用 Matlab 做字符串替换把数据文件名和阈值换掉。从那以后我每次算完 Liutex 都强制走一遍检查清单梯度计算用中心差分且边界单独处理、特征值容差按数据精度设、涡轴方向与涡量对齐、DAT 文件先小规模测试再批量写、Tecplot 宏先在单帧上验证再循环。这套流程跑顺了从速度场到涡结构出图基本可以做到半自动。希望帮到你。本文还有配套的精品资源点击获取