陳 槐 張 磊 王乃茹 朱立俊

* (南京水利科學研究院水文水資源與水利工程科學國家重點實驗室,南京 210029)

? (中國水利水電科學研究院水利部泥沙科學與北方河流治理重點實驗室,北京 100038)

引言

湍流兼具隨機性和相干性,其基本組成結構是各種尺度的渦旋,背景流場中存在大量隨機性的小尺度渦旋,還存在時間或空間上相關性很強的大尺度渦旋,渦旋被譽為流體運動的肌腱[1].由于渦旋與湍流、大氣現象、海洋物理、流體控制、升增減阻和氣動噪聲等方面關系密切,渦旋研究在航空航天、水利、海洋、動力機械和仿生學等領域工程應用中具有重要的意義[2].

渦旋被掩蓋在雜亂的湍動場中,需要利用特定的方法進行提取,最常用的是參數閾值法,其建立在局部速度梯度張量的基礎上,具有伽利略不變性,與參考系無關,主要分以下4 類[3]: (1)第二不變量Q＞0,表征了局部空間內的旋轉率大于應變率;(2)壓力Hessian 矩陣的第二特征值λ2＜ 0,表征了特定平面內的旋轉率大于應變率;(3)判別式Δ＞ 0,表征速度梯度矩陣存在復特征值;(4)旋轉強度λci＞ 0,為梯度矩陣復特征值的虛部,表征了渦旋旋轉的角速度.

雖然渦旋已經研究了幾十年,但鑒于三維渦旋結構的復雜性,已有文獻一般采用前述參數(Q,λ2,Δ及λci)某一固定閾值條件下的三維等值面體進行可視化,定性描繪三維渦旋結構的空間形態(tài)及其演化特征[4-17].當要定量描述三維渦旋時,為簡化難度,已有文獻主要是從平面二維流場中開展分析,量化渦旋的形態(tài)屬性,如數量、尺度及方向等(相關文獻較多,不再羅列,可參見文獻[3]引言).然而,這些渦旋二維形態(tài)的研究結果(僅列舉代表性文獻[18-24])很大比例是建立在Wu 等[24]的研究思路基礎上,選用不同壁面距離下的旋轉強度均方根對旋轉強度進行歸一化,使各展向平面內旋轉強度的概率密度函數曲線基本重合,由此得到渦旋識別的統(tǒng)一閾值.雖然這個方法看似閾值統(tǒng)一,實則閾值沿壁面距離發(fā)生改變,由其所得的結果必定與直觀認識的渦旋結構(可視化時采用固定閾值)存在較大的差異,還會在某種程度上“扭曲”讀者對三維渦旋結構的理解.

為了還原渦旋結構的真實三維形態(tài)屬性,本文直接從渦旋可視化時利用的三維等值面體出發(fā),研究三維渦旋的形態(tài)屬性,達到所見即所得的效果,而不必像先前那樣采用兩套不同流程(可視化時采用三維流場下固定閾值,定量計算時采用二維流場下變閾值).

本工作采用兩三角形相交快速檢測算法,提取渦旋多邊形,利用數量與形態(tài)參數(圓度、半徑、凸狀及縱橫比因子)刻畫三維渦旋的形態(tài)屬性,并研究渦旋等值面體的閾值與壁面距離對渦旋形態(tài)的影響規(guī)律.研究結果以期有助于進一步深入理解渦旋的形成、發(fā)展和衰亡以及渦旋之間、渦旋與其他因素(如泥沙顆粒)之間的相互作用規(guī)律和機理.

1 數據及方法

1.1 DNS 數據

本文采用Del-Alamo 等應用DNS 方法獲得的槽道湍流數據,該數據幾乎模擬出槽道湍流中所有的含能結構,被廣泛應用于相干結構的研究[5,22,25].槽道流的摩阻雷諾數為Reτ=uτh/ν=934,其中uτ為摩阻流速,h是槽道垂向半高,ν是動力黏滯系數,x,y及z分別對應流向(沿水流方向)、垂向(垂直于床面向上) 及展向(槽道展寬方向),詳細數值見表1,其中Δx和Δz是x和z方向的網格分辨率,Δyc是槽道中心的網格分辨率,Nx,Ny及Nz是對應的網格數;上標“+”表示用內尺度(uτ和ν)進行無量綱化,如Δx+=Δxuτ/ν.

表1 DNS 槽道湍流計算參數Table 1 Parameters of the DNS channel flow dataset

1.2 渦旋識別方法

在眾多的渦旋識別方法中,旋轉強度法應用廣泛,它能區(qū)分背景剪切,具有明確的物理意義(渦旋旋轉的角速度),可提取渦旋的旋轉平面及旋轉角速度,并在渦結構可視化時對閾值的要求較低[3].令Mij=?ui/?xj(i,j=1,2,3,其中1,2,3 分別表示沿流向x,垂向y,展向z方向)是三維不可壓流場中任意一點(x,y,z)的三維速度梯度張量,其特征方程為

其中Q=0.5tr(?MM),R=?det(M).利用卡丹公式求解此一元三次方程,當判別式Δ=(R/2)2+(Q/3)3＞0 時,方程存在一個實數根及一對共軛復數根,旋轉強度即是復數根的虛部,其表達式為

1.3 渦旋形態(tài)因子提取方法

利用移動立方體算法快速生成渦旋的三維結構等值面體,并設置三角形切面對其進行剖分,通過兩三角形相交快速檢測算法,提取渦旋三維結構等值面體的切線,構造封閉多邊形,由此計算三維渦旋的形態(tài)因子(如圓度、半徑、凸狀及縱橫比因子),量化渦旋三維結構形態(tài),具體步驟如下.

(1)在某一旋轉強度閾值下(一般采用某一高程平面或某一立方體內的歸一化旋轉強度或旋轉強度最大值衡量),利用移動立方體算法[26]快速生成渦旋的三維結構等值面體(以三角網格為單位,見圖1 中綠色等值面體).

圖1 渦旋等值面體與三角形切面Fig.1 3D vortex structure indicated by isosurface of swirling strength and a triangular cutting plane

(2)在某一高程y+處,設置xz面內的三角形切面對渦旋三維結構等值面進行切割(見圖1 中紅色切面),通過兩三角形相交快速檢測算法[27],提取大的三角形(切面)與每一小三角(等值面體的最小構成單元)間的相交線,由此得到非常多的切線段.

(3)利用渦旋的物理屬性(一個首尾相連的閉合多邊形),將步驟2 中的切線進行篩選,僅保留那些首尾各有一條其他切線相連的切線,并組合成封閉多邊形,從而提取出該三角形切面內的所有渦旋的位置信息(見圖2(a)),為避免誤差影響,僅保留那些多邊形的最小外接邊界框長、寬都大于網格最小尺寸(3.8+)的渦旋.

圖2 渦旋提取與形狀因子定義Fig.2 Vortex extraction and shape factors definition

(4)令步驟3 中提取的任意一個渦旋的節(jié)點坐標按逆時針方向排序為(xi,zi) (i=1,2,···,N),其中i為某個節(jié)點的序號,N為節(jié)點總數,則渦旋的形狀因子公式(見圖2(b))如下[28]

其中,O為圓度因子,面積濕周,r為半徑因子,c為凸狀因子,Ac是渦旋多邊形的凸包(Graham 算法)所圍成的面積,E為縱橫比因子,a與b分別是渦旋多邊形最小外接橢圓(Khachiyan 算法)的長軸與短軸.凸狀因子反映出封閉曲線的彎曲程度,如果曲線完全外凸(曲率大于0)則凸狀因子等于1,如果曲線有內凹(曲率小于0)則凸狀因子小于1.以67#渦旋為例,其圓度、半徑、凸狀及縱橫比因子分別為0.45,17.16,0.82 及0.32.

2 計算結果

2.1 結果驗證

圖3 給出了對數區(qū)y+=110 平面內本文結果與前人結果的對比.Gao 等[5]采用參數閾值法提取渦旋結構,計算步驟簡要概括為: 計算平面內所有網格點的旋轉強度值,標記出數值超過0.13 倍最大值的點,并篩選出局部極大值點作為渦旋核心,利用區(qū)域生長算法提取渦旋邊界,最終計算渦旋半徑.值得指出的是,當利用區(qū)域生長算法計算時,渦旋尺寸最少需要滿足9 個點(即一個渦旋中心點與8 個外圍點),所以Gao 等[5]提取出的渦旋最小半徑為,體現為概率密度函數(PDF)曲線在r+=9 時出現極大值.為與前人對比,本文也計算了T=閾值下的渦旋等值面體,并利用切面法提取渦旋半徑,由于會保留滿足渦旋多邊形最小外接邊界框長、寬都大于3.8 的渦旋,所以本文的渦旋半徑最小可以接近1.07).

圖3 本文結果與前人結果渦旋半徑PDF 對比圖(y+=110,樣本數量9 423 250)Fig.3 Comparison on PDF of vortex radius between this study and previous study (y+=110,the number of samples is 9 423 250)

由圖可知,由于采用區(qū)域生長算法,Gao 等[5]的計算結果無法捕捉到r+＜ 9 的渦旋,造成其渦旋半徑PDF 曲線沒有峰值,只在r+=9 處出現極大值;反觀本文,由于采用切線法,可以精確地捕捉出小尺度的渦旋,得到更符合物理意義的渦旋半徑PDF 曲線,渦旋半徑在r+=6.5 處出現峰值.此外,由于區(qū)域生長算法的限制及網格尺寸略大的原因,造成Gao 等[5]提取的渦旋半徑要略大于本文計算結果,但兩者的PDF 曲線在r+＞ 9 區(qū)間內變化規(guī)律基本一致.

2.2 閾值對形狀因子的影響

由于渦旋的三維結構等值面體是在某閾值下生成的,所以閾值的大小會顯著影響等值面體的大小.對于渦旋結構而言,旋轉強度閾值越大,等值面體所構成的渦管就會越細,相應的渦旋半徑就會越小.由于量化渦旋三維結構過于復雜,已有文獻中一般僅通過研究某一平面內的旋轉強度場來研究渦旋的屬性(如數量、尺度及方向等),閾值取為該平面內旋轉強度最大值的0.08～0.32 倍[5-6,12,16].為研究渦旋三維結構等值面體的閾值對渦旋形狀因子計算結果的影響,本文以對數區(qū)內y+=110 平面為例,研究閾值范圍條件下渦旋的數量與形狀因子的變化規(guī)律.

圖4 給出了閾值對渦旋提取數量的影響.將渦旋數量轉為無量綱化渦旋密度,即其中N是渦旋數量,x與z方向模擬區(qū)域尺寸分別為=3πuτ/ν.可知,閾值越大,渦旋提取密度越小.當閾值從0.05 增加至0.5,渦旋提取密度從10?2降至10?6,雖然閾值僅增加了一個量級,但渦旋密度卻降低了4 個量級.渦旋提取密度與閾值關系基本為對數關系,擬合曲線表達式為lgΠ=?8.25T?1.67,擬合相關系數R2=0.995 8.

圖4 閾值對渦旋提取密度的影響(y+=110)Fig.4 Influence of threshold on the number of extracted vortices(y+=110)

圖5 給出了閾值對渦旋形狀因子概率函數的影響,值得指出的是,圓度、半徑及縱橫比因子采用概率密度函數PDF 表示,而凸狀因子由于其PDF 曲線在各閾值下過于緊湊,只能采用概率累積分布函數CDF 將其區(qū)分.

圖5 閾值對渦旋形狀因子概率函數PDF 及CDF 的影響(y+=110)Fig.5 Influence of threshold on PDFs and CDFs of vortex shape factors (y+=110)

對于圓度因子,當閾值很小時,其PDF 曲線比較平坦,類似均勻分布,表明此時的渦旋主要以非圓形態(tài)存在;隨著閾值逐漸增大,PDF 曲線慢慢變得高聳(左偏分布),各閾值間的PDF 曲線也逐漸重合,表明渦旋形態(tài)逐漸穩(wěn)定,并主要以類圓形態(tài)存在(眾數為0.73).

對于半徑因子,當閾值較小時,等值面體接近渦旋邊緣,所以得到的渦旋半徑較大;當閾值較大時,等值面體接近渦旋核心,所以得到的半徑較小.PDF 曲線簇變化規(guī)律與圓度因子相似,隨著閾值增大,PDF 曲線逐漸高聳并重合(眾數在2 左右);兩者不同的是,半徑因子PDF 曲線是右偏分布.

對于縱橫比因子,其PDF 曲線隨閾值的增加并未顯著變化,僅右偏程度稍有減小,PDF 曲線簇基本重合在一起(眾數分布在[0.4,0.5]區(qū)間).可知,在各閾值下,渦旋多邊形的最小外接橢圓形狀統(tǒng)計規(guī)律基本一致,長、短軸長度比值在2～2.5 范圍內的橢圓最多.

對于凸狀因子,可以通過觀察CDF 曲線的斜率來反映其PDF.當閾值很小時,渦旋的凸狀因子數值主要分布在[0.4,1]區(qū)間內,且在[0.6,0.9]區(qū)間內PDF 基本呈均勻分布(CDF 曲線斜率相同),說明此時大部分渦旋多邊形都不飽滿,它們的輪廓存在凹陷,且凹陷程度都較大.隨著閾值的增加,渦旋凸狀因子數值分布區(qū)間逐漸變窄并向1 靠近,當閾值達到最大值時,其CDF 和PDF 曲線都類似脈沖函數,說明隨著閾值的增加,渦旋多邊形逐漸變得飽滿,輪廓幾乎不存在凹陷了.

圖6 給出了閾值對渦旋形狀因子均值的影響.隨著閾值的增加,圓度因子的均值先增加隨后基本在0.65 附近波動,表明渦旋形狀先不斷變得更圓,隨后基本穩(wěn)定不變;縱橫比因子的均值始終在0.5 上下波動,并未出現較大程度改變,表明渦旋外接橢圓長、短軸比值約維持為2;凸狀因子均值始終在增加并最終接近1,表明渦旋多邊形趨于飽滿;半徑因子均值遞減,當T＜ 0.2 時遞減速度較大,此后遞減速度基本維持不變.

圖6 閾值對渦旋形狀因子均值的影響(y+=110)Fig.6 Influence of threshold on mean values of vortex shape factors(y+=110)

2.3 壁面距離對形狀因子的影響

前人在可視化渦旋三維結構時采用參數閾值法(Q,λ2,Δ或λci),利用某一固定閾值下的等值面體定性描繪渦旋三維結構,而已有關于渦旋二維切面的形態(tài)屬性(數量、尺寸及方向)的定量研究結果都是建立在沿壁面距離變閾值的思路基礎上,這種變閾值提取的結果在某種程度上“扭曲”了我們對真實三維渦旋結構的理解,亟待研究固定閾值下的相關結果.

圖7 固定閾值()與各zx 平面內旋轉強度最大值()(y)的比值Fig.7 Ratio of fixed threshold (=0.04) to the maximum swirling strength in each zx-plane ()(y)

圖8 給出了壁面距離對渦旋提取數量的影響.隨著壁面距離的增加,渦旋密度在緩沖層內逐漸增加,并在對數區(qū)內達到最大值,此后沿程不斷衰減,極大值與極小值在全域相差約一個量級.與圖4 類似,圖8 中y+=[90,700]區(qū)間內渦旋密度與壁面距離也呈對數關系,lgΠ=?1.7×10?3y+?3.14,擬合相關系數R2=0.999 6,可能由于壁面和水面的影響,導致這兩個區(qū)域內關系不符合對數律.

圖8 壁面距離對渦旋提取密度的影響Fig.8 Influence of wall distance on the number of extracted vortices

圖8 還對比了本文定閾值結果與已有文獻變閾值結果,其中Chen 等[19]采用變閾值方法提取渦旋,將大于(平面內旋轉強度的均方根)的局部極大值點選出,且僅當其周圍8 個點同時大于時定義為渦旋.可見,兩條曲線差異較大,定閾值曲線快速遞減(極大值是極小值的10 倍以上),而變閾值曲線緩慢遞減(極大值僅是極小值的1.5 倍左右).此外,當y+＜ 580,定閾值渦旋密度大于變閾值,而當y+＞ 580,定閾值渦旋密度小于變閾值.原因主要有兩點: (1) 變閾值的在近壁區(qū)大于定閾值0.04,而在外區(qū)小于0.04;(2) Chen 等[19]用9 點法篩選渦旋(與2.1 節(jié)一致),而本文基于等值面體的提取方法捕捉渦旋,所以造成兩者的結果差異較大.

圖9 給出了壁面距離對渦旋形狀因子概率函數的影響.總的來說,圖9 與圖5 存在很多相似之處,都是隨著自變量(閾值或壁面距離)的增大,因變量(4 個參數)的PDF 或CDF 逐漸重合.渦旋形狀因子的概率函數在不同分區(qū)內表現的不同,以緩沖層內最為顯著.隨著壁面距離的增加,對于圓度因子,PDF 曲線從右偏變?yōu)樽笃植?眾數從0.4 變?yōu)?.9;對于半徑因子,PDF 曲線右偏分布,緩沖層內曲線較為平坦,眾數位于[15,30]范圍,對數區(qū)與尾流區(qū)曲線高聳且基本重合,眾數在8.5 附近;對于縱橫比因子,PDF 曲線在緩沖層內右偏分布,眾數約0.2,在對數區(qū)及尾流區(qū)內接近無偏分布,眾數約0.5;對于凸狀因子,從緩沖層到尾流區(qū),數值分布區(qū)間逐漸變窄并向1 靠近,各垂向平面的CDF 曲線與圖5 中相應比例閾值下的曲線基本一致.

圖9 壁面距離對渦旋形狀因子概率函數PDF 及CDF 的影響Fig.9 Influence of wall distance on PDFs and CDFs of vortex shape factors

圖10 給出了壁面距離對渦旋形狀因子均值的影響.隨著壁面距離的增加,圓度因子均值不斷增加(從0.43 增至0.74),但增速在緩沖層與對數區(qū)內較大,在尾流區(qū)內較小.

圖10 壁面距離對渦旋形狀因子均值的影響Fig.10 Influence of wall distance on mean values of vortex shape factors

與圓度因子類似,凸狀因子均值也隨壁面距離增加而增大(從0.84 增至0.96),增速也是先大后小,但數值上各區(qū)差異不大;縱橫比因子均值在y+=[5,70]區(qū)間內快速遞增,從0.35 增至0.54,此后基本維持不變;半徑因子在y+=[5,140]區(qū)間內快速遞減,從15 減至11.5,此后基本保持穩(wěn)定.

緩沖層內以準流向渦為主,傾角(渦管軸向矢量與xz面夾角)一般以10°～20°為主;對數區(qū)內發(fā)夾渦逐漸產生,并在尾流區(qū)內成為主要的渦結構,傾角主要以45°為主[3,30].由于緩沖層內渦旋結構傾角很小,相比其他兩個分區(qū),渦管切面的多邊形更狹長,所以緩沖層內計算的圓度、縱橫比因子數值都較小,而半徑因子數值較大.

圖10(b)給出了本文定閾值結果與已有文獻變閾值結果的對比,其中Herpin 等[22]和Chen 等[19]都采用變閾值方法提取渦旋,將大于(平面內旋轉強度的均方根)的局部極值選為渦旋中心,不同的是,Herpin 等[22]利用Oseen 渦模型擬合獲取渦旋半徑,而Chen 等[19]采用區(qū)域生長算法提取渦旋邊界并計算渦旋半徑.可以看出,Herpin 等[22]和Chen 等[19]由于采樣相同的變閾值,所以渦旋半徑都是隨壁面距離增大而增加的,差異僅是數值而已.本文采用定閾值得到渦旋半徑隨壁面距離先減小(與準流向渦角度有關)后保持不變.可見定閾值和變閾值得到的規(guī)律完全不同,已有的變閾值結果會讓讀者造成誤解,認為渦旋半徑隨壁面距離增大,其實如果閾值相同,渦旋半徑全域幾乎不變.

3 結論

本文基于渦旋等值面體研究了閾值與壁面距離對渦旋形狀因子的影響.主要通過旋轉強度提取渦旋信息場,利用移動立方體算法快速生成某閾值下的渦旋三維結構等值面體(以三角網格為單位),設置三角形切面,通過兩三角形相交快速檢測算法提取兩者的切線,結合渦旋物理屬性提取渦旋多邊形,由此計算渦旋形態(tài)因子(圓度、半徑、凸狀及縱橫比因子),并通過概率函數曲線與統(tǒng)計均值的變化,研究閾值與壁面距離對渦旋形狀因子的影響規(guī)律.值得指出的是,本文研究結論是在單一雷諾數條件下獲得的,其他雷諾數下的結果可能有一定的差異.本文得出主要結論如下.

(1)相比基于區(qū)域生長算法等的傳統(tǒng)渦旋提取方法,基于渦旋等值面體并利用切面法提取渦旋的方法,可以更精確地捕捉出接近網格尺度的渦旋,得出更符合物理意義的渦旋形態(tài)參數統(tǒng)計規(guī)律.

(2)量化了對數區(qū)y+=110 平面內閾值范圍條件下渦旋數量及形狀參數的變化規(guī)律.對于渦旋密度,閾值越大,渦旋密度越小,兩者呈對數關系 lgΠ=?8.25T?1.67.隨著閾值的增加,對于圓度因子,其PDF 曲線由平坦變得高聳,均值先增加隨后基本在0.65 附近波動,表明渦旋不斷變得更圓隨后保持類圓形態(tài);對于半徑因子,其PDF 曲線簇變化規(guī)律與圓度因子相似,但均值不斷減小;對于縱橫比因子,其PDF 曲線未顯著改變,均值在0.5 附近波動,表明渦旋外接橢圓長、短軸比值約維持為2;對于凸狀因子,數值分布區(qū)間逐漸變窄并向1(類似脈沖函數)靠近,說明渦旋逐漸變得飽滿.

(4)闡明了全域定閾值與變閾值條件下渦旋數量與半徑變化規(guī)律的差異.對于渦旋密度,定閾值與變閾值結果都隨壁面距離先增大后減小,但變閾值導致渦旋密度在尾流區(qū)衰減不符合對數律.對于變閾值,已有文獻中渦旋半徑隨壁面距離增大而增加;對于定閾值,渦旋半徑隨壁面距離快速減小(與流向渦有關)后基本保持不變.這種人為刻意的沿壁面距離變閾值的渦旋屬性提取的結果,“扭曲”了我們對真實三維渦旋結構的理解,誤認為在遠離壁面的過程中渦核會變大.其實只要閾值一定,渦核尺寸就維持不變,遠離壁面的過程中只是高能渦旋不斷破碎并將傳遞能量給低能渦旋,減小的只是高能渦旋的數量.