簡介該資源為Liutex渦定義與第三代渦識別方法的入門與實現(xiàn)資料包面向流體力學(xué)、渦動力學(xué)方向的研究人員及工程技術(shù)人員幫助讀者理解Liutex理論并快速上手Matlab與Tecplot的可視化實現(xiàn)。壓縮包共6個文件包括3個Matlab腳本如getR.m、writeLiutex.m等用于計算與輸出Liutex場、1份PDF文獻介紹、1份PPTX講解材料及1個公式文件整體大小6.67MB內(nèi)容精簡但覆蓋從原理到編碼的關(guān)鍵環(huán)節(jié)。已有3582人瀏覽學(xué)習(xí)較為適合具有一定流體力學(xué)基礎(chǔ)、希望掌握新代渦識別方法的讀者。資料由作者針對網(wǎng)上盜賣行為重新免費上傳無需積分即可下載。通過閱讀PDF與PPT可系統(tǒng)了解Liutex定義與優(yōu)勢借助Matlab腳本可對照實現(xiàn)核心算法并進一步結(jié)合Tecplot完成結(jié)果可視化從而在科研或工程項目中直接應(yīng)用這一較前沿的渦識別工具。1. Liutex 不是玄學(xué)為什么我用 Matlab 算、用 Tecplot 看做 CFD 后處理的人大多會遇到同一種尷尬算完流場用 Q 準(zhǔn)則或 λ2 準(zhǔn)則挑渦結(jié)構(gòu)結(jié)果射流剪切層、壁面邊界層、甚至純剪切縫隙里都飄著一大片“渦”可這些區(qū)域根本沒有旋轉(zhuǎn)核心云圖出來后完全沒法跟實驗流動顯示對上。Liutex 方法就是來補這個空子的它把速度梯度張量拆成剛性旋轉(zhuǎn)和變形兩類貢獻只保留真正的剛性旋轉(zhuǎn)強度 R再輸出矢量。這條路我用 Matlab 做逐點計算、Tecplot 做后處理出圖腳本全流程跑通不依賴商業(yè) CFD 求解器內(nèi)部函數(shù)。適合手上只有速度場結(jié)果、想快速驗證渦結(jié)構(gòu)的仿真工程師和研究生。2. Liutex 的數(shù)學(xué)原理與計算路線從速度梯度到旋轉(zhuǎn)軸2.1 為什么不再依賴 Q 準(zhǔn)則和 λ2 判據(jù)傳統(tǒng)判據(jù)的通病在于它們都建立在“渦量大于變形就等于有渦”這個近似上?？疾焖俣忍荻葟埩?u的對稱部分 S 和反對稱部分 ΩQ 準(zhǔn)則計算0.5*(||Ω||2 - ||S||2)λ2 準(zhǔn)則看的是 S2Ω2 的第二特征值。這兩個判據(jù)在純剪切層里的響應(yīng)非常糟糕純剪切場的渦量和變形量可以同時很大Q 值照樣偏高于是把一大片剪切區(qū)誤報成渦。我最早用 Q 準(zhǔn)則時判斷“哪里是渦”基本靠猜后來才意識到這是剪切污染。Liutex 的計算思路完全不同。它不比較渦量與變形的相對大小而是先問一個問題這個局部區(qū)域是否存在一條旋轉(zhuǎn)軸能讓流體質(zhì)點繞它做剛性旋轉(zhuǎn)判斷依據(jù)是速度梯度張量的特征值結(jié)構(gòu)。只要存在一對共軛復(fù)特征值和一條唯一的實特征值對應(yīng)的旋轉(zhuǎn)軸再滿足強旋轉(zhuǎn)判別式就認為該點是渦。這個門檻把剪切、拉伸這類非旋轉(zhuǎn)變形濾掉了所以它輸出的等值面通常比 Q 準(zhǔn)則小兩三圈正好對應(yīng)流體繞某條軸轉(zhuǎn)動的核心。2.2 特征值分解與 Liutex 強度的計算完整算法分五步。第一步計算每一點的速度梯度張量?u第二步求?u的三個特征值判斷是否存在一個實特征值和一對共軛復(fù)特征值如果三個特征值都是實數(shù)該點直接置零第三步取實特征值對應(yīng)的特征向量單位化后作為旋轉(zhuǎn)軸方向 r第四步計算渦量 ω 在 r 上的投影ω_r ω·r并取復(fù)特征值的實部 λ_cr第五步求剛性旋轉(zhuǎn)強度 R。公式是R sign(ω_r) * sqrt( max(ω_r2 - 4*λ_cr2, 0) )Liutex 向量就是R·r。公式里的4*λ_cr2項就是把速度梯度張量中由拉伸/剪切帶來的變形貢獻扣除剩下才是繞 r 軸的剛性旋轉(zhuǎn)分量。這套計算在 Matlab 里可以做成一個獨立函數(shù)方便后續(xù)逐點調(diào)用。function [Rvec, Rmag] liutex_grad(D, omega) % D : 3x3 速度梯度張量 % omega : 1x3 渦量向量 % 返回 Liutex 向量 Rvec 與標(biāo)量強度 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; % 三個實根或退化情況不產(chǎn)生剛性旋轉(zhuǎn) 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; % 統(tǒng)一旋轉(zhuǎn)軸方向避免相鄰點符號跳變 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邏輯說明每個網(wǎng)格點單獨調(diào)用這個函數(shù)返回零向量說明該點不是渦。n_real用容差判斷而不是imag(lam) 0因為浮點計算里純實數(shù)根的虛部也可能出現(xiàn) 1e-16 量級的殘留嚴格相等會把好點全部漏掉。旋轉(zhuǎn)軸方向統(tǒng)一那一步不是可選項特征向量只有方向沒有正負同一物理渦在相鄰網(wǎng)格點可能一個取正一個取負不做對齊則輸出 R 向量會連續(xù)倒向。參數(shù)說明tol取張量 F 范數(shù)的 1e-8 倍比固定閾值如 1e-10 更抗量綱影響壓氣機流場速度量級大、梯度張量數(shù)值大固定閾值很容易把所有特征值都判成實根lam_cr取第一個虛部超過容差的特征值的實部一對共軛復(fù)特征值實部相同取哪一個結(jié)果一樣。R 和 λ_cr 的單位都是 1/s與渦量一致所以后處理時可以直接和渦量場對比量級。2.3 旋轉(zhuǎn)軸方向統(tǒng)一關(guān)鍵數(shù)值細節(jié)特征向量作為旋轉(zhuǎn)軸方向時存在一個天然歧義r和-r是同一個方向。如果不在代碼里做統(tǒng)一輸出 R 向量會在相鄰網(wǎng)格點之間出現(xiàn)隨機反向。我一般用dot(r, omega) 0來判斷讓 r 與局部渦量方向的點積非負。這樣旋轉(zhuǎn)軸方向在空間上連續(xù)R 向量在流場中的指向才有物理意義。極少數(shù)退化點上渦量接近零或 r 與 ω 接近垂直點積接近零符號判定會抖動。這類點本身 R 也趨于零對渦核等值面影響可以忽略。不要在代碼里把這類點特殊處理成 NaNTecplot 顯示時會留下空洞反而干擾后續(xù)處理。3. Matlab 實現(xiàn)從 CFD 速度場到 Liutex 矢量的完整流程3.1 讀入結(jié)構(gòu)化網(wǎng)格速度數(shù)據(jù)假設(shè)已經(jīng)有計算出的結(jié)構(gòu)化網(wǎng)格速度場無論是自己寫的求解器還是從商業(yè)軟件導(dǎo)出的 Tecplot ASCII 文件這一步都先把數(shù)據(jù)讀到 Matlab 四維數(shù)組中。Tecplot ASCII 的 FPOINT 格式比較簡單數(shù)據(jù)按 I 方向最快、J 方向其次、K 方向最慢排列每個點按順序?qū)?X、Y、Z、U、V、W。我習(xí)慣先掃描到 ZONE 關(guān)鍵字再開始讀數(shù)避免不同軟件導(dǎo)出的表頭行數(shù)不一致。function [x,y,z,u,v,w] read_tecplot_ascii(fname, ni, nj, nk) % 讀取 Tecplot ASCII 點序數(shù)據(jù) % 變量順序固定為 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 的行找到后文件指針已經(jīng)停在數(shù)據(jù)起點textscan從下一行開始讀取所有浮點數(shù)。reshape成 6 列是為了把每個點的分量拆開再用reshape還原成三維數(shù)組時第一維對應(yīng) I 方向和 Tecplot 的 FPOINT 寫入順序一致。參數(shù)說明ni、nj、nk必須與 Tecplot ZONE 行里的 I、J、K 一致否則 reshape 會報錯或靜默錯位。讀取后建議先對 x、y、z 做一次單調(diào)性檢查diff(x(:,1,1))出現(xiàn)變號說明數(shù)據(jù)順序理解反了。3.2 核心函數(shù)逐點計算 Liutex 向量讀入速度分量后用 Matlab 自帶的gradient函數(shù)計算九個梯度分量再組裝成速度梯度張量 D調(diào)用上一章寫好的liutex_grad。注意gradient的三參數(shù)版本要求網(wǎng)格等間距結(jié)構(gòu)網(wǎng)格且 dx、dy、dz 為常數(shù)時可以直接用。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 的三個方向?qū)?shù)組裝順序要和liutex_grad里的eig(D)對應(yīng)。渦量用旋度公式直接算符號約定是右手系如果數(shù)據(jù)來自左手坐標(biāo)系R 會整體反向需要在讀入階段先檢查坐標(biāo)系。參數(shù)說明gradient在邊界上用單側(cè)差分邊界層第一層網(wǎng)格的梯度和 R 會偏大這是后處理時最需要注意的誤差來源后面避坑章節(jié)會展開。三重循環(huán)在五百萬網(wǎng)格點上大概要跑幾分鐘到十幾分鐘主要開銷在eig如果網(wǎng)格過億可以把外層循環(huán)改成parfor但liutex_grad內(nèi)部不要共享變量。3.3 導(dǎo)出 Tecplot 可直接打開的數(shù)據(jù)文件計算完成后把速度原始場和 Liutex 場一起寫入新的 Tecplot ASCII 文件。這樣導(dǎo)入 Tecplot 后既可以同時顯示原始速度云圖和 Liutex 渦核又能在同一份數(shù)據(jù)里做截面和向量圖不用來回切換文件。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 TLiutex, I%d, J%d, K%d, FPOINT\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 的 FPOINT 順序一致因為數(shù)據(jù)寫入時 I 方向最快。變量列表里 Rx、Ry、Rz 是 Liutex 向量三個分量Rmag 是標(biāo)量強度Tecplot 里既能用 Rmag 做等值面也能用 Rx/Ry/Rz 畫旋轉(zhuǎn)軸向量。參數(shù)說明%.6e足夠后處理繪圖精度不建議用%f因為速度量級跨度大固定小數(shù)位會丟失小量信息。老版本 Tecplot 打開時報 zone 類型錯誤可以在 ZONE 行加ZONETYPEORDERED。寫完文件后用load(liutex_field.dat)讀回來檢查第一列是否單調(diào)遞增這是最廉價的驗證方式。4. Liutex 實現(xiàn)的常見坑現(xiàn)象、原因與排查4.1 旋轉(zhuǎn)軸符號跳變導(dǎo)致 R 向量方向混亂現(xiàn)象輸出的 R 向量在相鄰網(wǎng)格點之間突然完全反向Tecplot 里向量圖呈鋸齒狀等值面邊界出現(xiàn)隨機毛刺。原因特征值分解給出的特征向量方向不確定r和-r是同一根軸逐點計算時不同網(wǎng)格點可能落到不同符號分支上。解決在liutex_grad里強制dot(r, omega) 0時取反。如果做完這步還有個別點跳變檢查這些點的渦量是否接近零接近零的區(qū)域 R 也接近零不影響渦核判斷不用處理。4.2 特征值虛部判斷失誤把剪切錯認成旋渦現(xiàn)象層流庫埃特流這種純剪切場Rmag 仍然和真實渦核同量級等值面密密麻麻。原因代碼里用了imag(lam) 0判斷實根浮點誤差讓純實數(shù)根的虛部出現(xiàn) 1e-16 量級殘留被誤判成復(fù)數(shù)特征值于是走進了旋轉(zhuǎn)分支。解決用相對容差tol 1e-8 * norm(D, fro)并且只要n_real ~ 1就返回零向量。更保守的做法是同時要求disc 0才給非零 R根號內(nèi)非正的區(qū)域一律按非渦處理。4.3 sqrt 內(nèi)判別式出現(xiàn)負值現(xiàn)象個別網(wǎng)格點 Rmag 突然出現(xiàn) NaNTecplot 等值面在這些點附近斷開。原因ω_r2 - 4*λ_cr2理論上是非負的但速度場經(jīng)過插值或數(shù)據(jù)格式轉(zhuǎn)換后局部梯度會引入噪聲導(dǎo)致根號內(nèi)出現(xiàn) -1e-10 量級的負數(shù)另一種可能是該點本身是強變形主導(dǎo)兩個共軛復(fù)特征值的實部過大不滿足強旋轉(zhuǎn)條件。解決用max(disc, 0)保護同時在腳本里統(tǒng)計disc 0的點數(shù)。如果負值點占比小于總點數(shù)的 0.1%可以直接忽略如果超過 1%說明輸入速度場沒有經(jīng)過良好插值或者網(wǎng)格變換時梯度算錯了不要靠改公式掩蓋。4.4 Tecplot 打開數(shù)據(jù)錯位變量數(shù)與 FPOINT 不匹配現(xiàn)象導(dǎo)入 liutex_field.dat 后坐標(biāo)云圖出現(xiàn)條紋或者某個變量顯示出來的數(shù)值明顯不對。原因fprintf 的占位符數(shù)量與 VARIABLES 聲明的變量數(shù)不一致常見的是少寫一個分量或多寫一個分量Tecplot 會按固定列寬讀取錯位后整個場都亂掉。解決檢查 VARIABLES 字符串里有幾個變量名fprintf 里每行就寫幾個%e。我一般寫完后立刻用textscan讀回來對比 Rmag 總量是否與計算時一致。另一個隱蔽問題是文件被 Excel 打開后另存換行符和空格都會被改Tecplot 可能讀不出正確的點數(shù)所以別讓編輯器碰 .dat 文件。4.5 邊界差分誤差讓 R 值虛高現(xiàn)象壁面附近出現(xiàn)一層條帶狀偽渦Rmag 在貼近邊界的第一層網(wǎng)格上突然增大向內(nèi)卻迅速下降。原因gradient在邊界上只能用單側(cè)差分一階精度比內(nèi)部中心差分低一到兩個數(shù)量級速度梯度過大導(dǎo)致特征值誤判R 在邊界被顯著放大。解決先分清是物理邊界還是人為截斷邊界。物理壁面處建議把計算域收縮一個網(wǎng)格跳過最外層兩點再統(tǒng)計數(shù)據(jù)量足夠時也可以在計算前對速度場沿壁面法向做一次五點平滑但要注意不要過度平滑導(dǎo)致真實渦核消失。邊界層的 R 虛高是所有特征值類渦判據(jù)的共性問題Q 準(zhǔn)則同樣存在不是 Liutex 特有的。5. 驗證腳本正確性的三個解析場與 Tecplot 出圖技巧5.1 用解析速度場做量級驗證拿到新數(shù)據(jù)直接出圖之前花十分鐘用三個解析速度場把腳本跑一遍能避免大量后期返工。表格是三個驗證用例用例速度場定義預(yù)期 Liutex 標(biāo)量 R剛體旋轉(zhuǎn)u-Ωy, vΩxR2Ω純剪切uεy, v0R0剛體旋轉(zhuǎn)拉伸u-Ωyεx, vΩx-εyR2√(Ω2-ε2)剛體旋轉(zhuǎn)用來驗證主量級是否正確純剪切用來驗證是否會被誤判成渦旋轉(zhuǎn)加拉伸用來驗證變形貢獻是否被正確剝離。第三組要求 |ε| |Ω|否則沒有剛性旋轉(zhuǎn)主導(dǎo)區(qū)R 應(yīng)當(dāng)返回零。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; % 后續(xù)復(fù)用 3.2 節(jié)的梯度與循環(huán)代碼 % 取中心點 (nx/2, ny/2) 的 Rmag 與 2*sqrt(Omega^2-Epsilon^2) 對比邏輯說明網(wǎng)格中心點遠離邊界差分誤差最小適合做定量對比期望值與數(shù)值解誤差應(yīng)小于 1e-8。如果誤差很大優(yōu)先檢查組裝 D 時梯度分量放的行列位置是否反了這個錯誤會把旋轉(zhuǎn)方向算錯但不影響特征值實部非常隱蔽。參數(shù)說明Ω 取 1.0 是為了讓量級便于檢查剛體旋轉(zhuǎn)的渦量是 2Ω所以 Rmag 期望值是 2.0不需要額外換算。5.2 Tecplot 里畫渦核等值面和旋轉(zhuǎn)軸向量數(shù)據(jù)導(dǎo)入 Tecplot 后我一般先做一個 Rmag 的等值面閾值取全場最大值的一半快速看渦核空間位置再用 Rx、Ry、Rz 畫 Vector確認旋轉(zhuǎn)軸方向與渦量方向大致一致。等值面不要直接拖 Contour 默認值先看 Rmag 分布直方圖找到有明顯拐點的閾值避免被邊界偽渦拉高最大值后等值面什么都顯示不出來。一個實用的步驟是在 Tecplot 里創(chuàng)建兩個切片一個垂直于旋轉(zhuǎn)軸看截面渦核一個沿旋轉(zhuǎn)軸看延展長度。如果渦核在切片上呈對稱圓形說明旋轉(zhuǎn)軸估計正確如果被拉成橢圓形說明該區(qū)域存在較強拉伸需要回到 Liutex 向量而不是只盯著 Rmag。有一次我圖省事拿到新的模擬數(shù)據(jù)直接跑腳本出等值面結(jié)果壁面全是偽渦結(jié)構(gòu)。排查了半天才發(fā)現(xiàn)是特征值判斷沒加容差把剪切區(qū)的實數(shù)根誤判成了復(fù)共軛根。從那以后我每次換網(wǎng)格格式、換求解器都強制先用這三個解析場把腳本過一遍再上真實數(shù)據(jù)。希望幫到你。本文還有配套的精品資源點擊獲取