渦識別技術及其在船用螺旋槳伴流場中的應用1)

楊 琳 鄭 興

(哈爾濱工程大學船舶工程學院,哈爾濱 150001)

引言

旋渦的概念與流體動力學的主題一樣古老,但是仍然缺乏公認的旋渦定義.關于湍流中旋渦是由什么構成的,長期以來存在相當大的困惑.通常可以將湍流中的擬序結構(coherent structure,CS)視為旋渦,實踐中必須首先確定湍流中動態顯著的大尺度旋渦區域,而這又需要對旋渦進行客觀定義.一般來說旋渦的客觀定義應允許使用渦動力學概念來推導和解釋CS 的形成與演化,探索CS 在湍流現象中的作用,為湍流現象開發可行的湍流模型和控制策略.湍流被視為渦流細絲的纏結,并且使用渦動力學的概念可以很好地解釋湍流物理學的許多方面[1].通過描述CS 的形成、演化及其與背景湍流的耦合過程,有望了解各種湍流現象,例如夾帶(entrainment)和混合、傳熱和傳質、化學反應和燃燒、阻力以及空氣動力學噪聲的產生,而且還可以對湍流進行可行的建模[2].

直覺上渦流被認為是一個管,其表面由渦線組成[3],然而渦流管的存在并不意味著渦流的存在,例如層流管流中的渦流管在任何情況下都不是渦流.Lugt[4]將旋渦定義為圍繞一個公共中心旋轉的大量物質顆粒;Chong 等[5]將其定義為 ?u的復特征值區域;Hunt 等[6]將其定義為一個區域,同時包含 ?u的正第二不變量和低壓;Jeong 等[7]根據對稱張量的特征值定義它,其中S和 ? 分別是速度梯度張量 ?u的對稱和反對稱部分.

研究已經發現湍流剪切流主要由空間連貫的和隨時間變化的旋渦運動所主導,通常表現為CS[8-10].在物理實驗和數值模擬中已經提出了幾種用于CS 生成的條件采樣技術[11-16].在過渡流中,CS 生成相對簡單,因為CS 的發生在時間和空間上具有規律性,諸如速度之類的參考信號可以用作生成的觸發信號[17-18];但是在充分發展的湍流中,例如射流(jets)、伴流(wakes)和混合層(mixing layers)的遠場區域,要得到動力學意義的旋渦結構,瞬時渦度場是必不可少的[19-20],遺憾的是即使是瞬時渦度場也不足以揭示湍流邊界層中的CS[21].

認識到黏性流體中旋渦的大小取決于所選識別器的閾值,并且大量研究表明湍流中渦旋 C S 的核在空間中定位良好,故選取僅識別渦旋核.渦旋核的存在一般要滿足兩個要求: (1)渦旋核必須具有凈渦度及凈環量;(2)識別出的渦旋核的幾何形狀應該是伽利略不變的.遺憾的是這些要求不會導致單一的識別方法,僅將這些要求用作對可能潛在識別方法的初步檢查.

國內對渦識別技術的理論研究及應用相對較少,在百度學術上搜到的近年公開發表的關于這方面的研究見文獻[22-24];楊越[22]總結了湍流中目前主要的渦識別方法,提出了渦識別和相關渦動力學的重要科學問題,并展望了湍流中渦動力學與相關湍流預測模型的發展趨勢;王義乾等[23]講到了第三代渦識別方法及其應用綜述;劉超群[24]回顧了渦定義和渦識別方法的發展歷史,著重介紹了美國德州大學阿靈頓分校團隊及其合作者在渦科學和湍流研究的一些最新學術創新成果.國外最新公開發表的關于渦識別方面的研究見文獻[25-35].盡管目前存在一些渦識別技術的研究,也有人提出一些新的識別方法,但是在現實中大部分人還是沒有很好地去理解什么是渦及如何去識別它們這兩個問題;且渦識別技術里包含許多深層次的數學及物理認識,還遠沒有到下結論的時候,值得更多的研究者去深入探討.

本研究旨在船用螺旋槳的伴流場中,從理論和實踐上尋找一種可行且有效的渦識別方法.如何識別一船用螺旋槳伴流場中的旋渦？本文作者圍繞這一問題進行了深入的探討,從理論上研究了六種渦識別方法,如表1 所示,并討論了各種識別方法的優缺點.為了更好地說明渦識別技術的合理性,本文還使用了一些解析解、Burgers 渦流解和Lamb-Oseen渦流解.

表1 渦識別方法Table 1 Methods of vortex identification

1 研究方法

1.1 控制方程與數值方法

伴流場基于納維?斯托克斯(N ?S)方程組的數值模擬,采用一個局部動態k方程大渦模型(large vortex simulation,LES)[36],在物理空間中,笛卡爾坐標系x=(x1,x2,x3),下速度場為u=(u1,u2,u3),質量守恒及動量守恒方程(不計體力)分別為

其中,p為靜壓,ν 為運動黏度并假設其為常數.根據不同的時空分辨率和動力學描述,LES 僅用于直接計算空間中的低頻模式,如果忽略濾波器和空間導數之間的換相誤差,過濾后的質量守恒及動量守恒方程分別為[37]

其中亞格子尺度(subgrid scale,SGS)應力張量 τ 和亞格子尺度動能kSGS方程分別為

其中,νSGS為亞格子渦黏度,ε 為亞格子動能耗散率.

LES 基于大尺度和小尺度之間的分離,其中直接模擬具有較高能量的大尺度湍流結構(各向異性湍流),而小尺度湍流則使用亞格子尺度模型進行建模.小尺度湍流波動本身被平滑,并通過渦流黏度和擴散率假設進行建模.特征尺寸大于截止長度的尺度稱為大尺度或分辨尺度,其他尺度稱為小尺度或亞格子尺度.SGS 通過統計亞格子尺度模型將其包括在內,SGS 的精確建模是可能的,因為在較小的尺度(低于過濾器寬度)下,湍流可以被認為是各向同性的,與流動類型和邊界條件無關.LES 就是對過濾后的 N ?S 方程中非線性項進行分解,然后對未知項進行建模.

本文研究了一種五葉右旋船用螺旋槳DTMB 4381[38](變槳距、無側斜、無尾斜)的伴流場.在無量綱半徑r′=r/R=0.2 ～1 處,葉片在與螺旋槳同軸的圓柱面上的截面為翼型,如圖1 所示,其中R是螺旋槳盤半徑,槳轂呈錐形.相對弦長c′(用螺旋槳直徑D無量綱化,即c′=c/D)和扭角 θT的分布如圖2 所示.為了消除數值水筒的阻塞效應,計算域的設置如圖3 所示,圖中也給出了流域的初始邊界條件.為了能捕獲所關心的伴流渦,網格劃分示意如圖4 所示,計算域內不同區域網格的大小見表2.注意,后面“結果與分析”中給出的螺旋槳伴流場的渦識別結果是在此螺旋槳的設計工況下得到的,即計算域入口處的來流速度U=4.5 m/s,螺旋槳的進速系數J=0.9,對應槳的轉速約1500 r/min.

圖1 螺旋槳幾何Fig.1 Geometry of blade,hub and shaft

圖2 葉片沿徑向相對弦長和扭角分布Fig.2 Distributions of relative chord length and twist angle along radial direction

圖3 計算域及初始邊界條件Fig.3 Computational domain and initial boundary conditions

圖4 網格劃分示意Fig.4 View of meshing

表2 計算區不同域網格尺度Table 2 Overall cell scales in different regions of computational domain

1.2 渦識別方法

在笛卡爾直角坐標系x=(x1,x2,x3)=(x,y,z) 下,若速度向量為u=(u1,u2,u3)=(u,v,w),則速度梯度張量T為

速度梯度張量T的對稱部分S和反對稱部分 ? 分別為

對稱張量S稱為應變率張量,其主對角線上的分量代表線應變變化率,其余分量代表角變形率;線應變率反映拉壓,角變形率反映剪切,兩者相互關聯;單位體積流體微團在單位時間內的體積變化率為張量S的第一不變量,流體微元體的體積只受線應變率的影響,而不受角變形率影響,角變形率僅引起剪切或角度變形.反對稱張量 ? 中僅包含三個不同的分量,對應一個矢量 ζ=rotu/2,描寫微團的瞬時角速度;由于速度的旋度與微團瞬時角速度之間的這種聯系,它可以作為流體渦旋運動強度的度量,稱之為渦強度或渦度,即有 ω=rotu=2ζ.

1.2.1 局部最小壓力標準

Robinson[21]使用DNS 數據庫顯示: 低壓標準有效地捕獲了湍流邊界層中的旋渦結構,但是通常無法指定適當的壓力水平來識別流動中的所有旋渦區域.最小壓力標準的基礎是: 在受離心力和壓力平衡(所謂的回旋平衡)的情況下,壓力趨向于在回旋運動的軸上具有局部最小值,此標準僅在定常無黏平面流中嚴格正確.壓力可能在某個點的所有方向上都具有最小值(所謂在三維度上局部最小壓力的概念),或者僅在垂直于渦旋軸的平面(例如Burgers渦)中具有最小值,僅討論平面中的最小值,此時對應的限制較少,且可能也包括前一種情況.最小局部壓力標準的不足之處可以通過四個簡單例子進行說明: (1)在非定常無旋流動中,存在一個明確定義的壓力最小值,考慮具有駐點的非定常無旋軸對稱流動

積分 E uler 方程,可得壓力為

控制的(此處的下標逗號表示偏微分,并且對指標使用求和約定),因此壓力固有地具有比渦度更大的尺度.

1.2.2 流線或跡線標準

Lugt[4]提出使用封閉或螺旋路徑來檢測旋渦,但是流線或跡線顯然不能滿足伽利略不變性.該檢測的另一個嚴重不足之處是,當旋渦由于諸如配對(pairing)、撕裂(tearing)等非線性過程分解,粒子在渦流的生命周期內可能無法完成圍繞旋渦中心的完整旋轉,沒有閉合的跡線.在渦流相互作用的重新連接區域,特別是在高R e時,由于快速改變的旋渦線拓撲結構,流線和跡線即使在最佳的參考系中也可能會高度扭曲,從而逃避被識別,而這些被忽略的區域可能在動力學上反而變得更加重要[41].

1.2.3 渦度值標準

渦度大小(| ω|)已被廣泛用于推導CS 并代表渦旋核[42-43],這種方法在自由剪切流中相當成功,但可能并不總是令人滿意,因為如果背景切變與渦旋內的渦度大小相當的話,| ω| 不能在剪切流中識別出渦旋核.Lugt[4]研究表明在平面邊界壁流中 | ω| 最大值和最小值僅出現在壁面上,Jimenez等[44]研究表明在湍流邊界層中最大值 |ω| 出現在高速條帶區域附近,然而在這兩種情況下,邊界附近的流動都以剪切為特征,沒有表現出旋動,故 | ω| 也不是邊界層中渦識別的合適標準.即使在自由剪切流中,從 | ω| 表面識別旋渦也存在潛在的困難,例如,渦片(vorticity sheet)不是渦旋(vortex),即使它可能具有很大的渦度值.

1.2.4 張量T的復特征值

由于速度梯度張量 ?u是伽利略不變量,Chong等[5]使用速度梯度張量 ?u的特征值對流場中任何點周圍的局部流線模式進行分類(參考系隨該點的速度移動),并提出渦核是一個具有復特征值的區域,復特征值表示局部流線樣式在隨點移動的參考系中是閉合的或螺旋形的.?u的特征值滿足特征方程(10),其中P=ui,i=0 (對不可壓縮流)、和它們分別是速度梯度張量的三個不變量,當判別式

為正值時,即存在復數特征值.

1.2.5 張量T的第二不變量Q或運動渦度數Nk

Hunt 等[6]定義渦流是具有 ?u的正第二不變量(Q>0)的區域,其附加條件是壓力低于環境值.Q定義為

表示剪切應變率和渦度大小之間的局部平衡,可以很容易地看出,與 | ω| 標 準不同,Q在壁面上消失了,即在壁面上的剪切應變率和渦度具有相同的大小.一般壁面上的速度梯度可以寫為

使得在壁面上有Q=0,這必然意味著 ?=0,因為在壁面處R=0,故基于 ? 和Q的定義不會出現 | ω| 定義無法正確地表示壁面附近的渦旋運動問題.Truesdell[45]引入了運動渦度數N k來測量旋渦質量,而不是 ∥?∥給出的局部旋轉速率,他將N k定義為

式中,N k是 | ω| 的逐點度量,通過應變率的范數進行歸一化處理,從而給出了旋渦質量,它與旋渦強度無關,例如N k=∞ 和N k=0 分別對應于固體旋轉和無旋運動,與 | ω| 的值無關.Melander 等[46]在研究渦核動力學時,將軸對稱旋渦圓柱的核心確定為N k>1 的最大連通空間區域.從方程(18)可知,N k>1 的區域與Q>0 的區域相同,但是由于N k通過應變率的大小進行了無量綱化,因此N k峰值忽略了旋渦動力學意義,換句話說只要旋渦質量相同,N k就不會區分不同渦度大小或環量的渦.Q也可以解釋為 N ?S 方程中的壓力源項,根據壓力泊松方程(?2p=2ρQ)和其最大原理,如果Q>0,最大壓力僅出現在邊界上;如果Q<0,最小壓力僅出現在邊界上.但是Q>0,并不一定意味著在該區域內會出現壓力最小值,最小壓力可以出現在Q>0 區域的邊界上,故在Q>0 的區域和最小壓力的區域之間沒有明確的聯系.從這個意義上說,Hunt 等[6]和Melander 等[46]使用的定義嚴格來說并不完全相同,即使在大多數情況下它們是等效的.從方程(15) 可以看出,條件Q>0 比?>0更具限制性,盡管哪個定義更合適尚無先驗.

1.2.6 張量S2+?2的特征值 λ2

如前所述,盡管不能將壓力局部最小值標準用作渦核識別的一般檢測標準,但是它卻為尋求合理的渦核識別提供了起點.局部最小壓力的存在并不表示渦核的存在,導致它們之間這種不一致性的主要原因有兩個: (1)非定常應變,它可以在不涉及渦旋運動的情況下產生最小的壓力;(2)黏性作用,它可以消除渦流運動中的最小壓力.通過簡單地丟棄這些影響,可期望獲得關于旋渦存在的更好標準.由于有關局部壓力極值的信息包含在壓力的Hessian(p,ij)中,考慮p,ij的方程,取 N ?S 方程的梯度可得

其中,a i,j為加速度梯度,p,ij對稱,a i,j可以被分解為對稱和反對稱部分

反對稱部分是渦度輸運方程,對稱部分可表示為

在平面中,局部最小壓力的出現需要張量p,ij的兩個正特征值.式(21)左側由于第一項表示非定常無旋應變,第二項代表黏性作用,將不考慮,僅通過考慮來確定是否存在由于旋渦運動引起的局部壓力最小值,并認定渦核為具有的兩個負特征值的連接區域.由于是對稱的,因此它僅具有實特征值,即有 λ1≥λ2≥λ3,則可定義渦核內 λ2<0 的標準,可得[45]

表3 基于負 λ2 和正 Q 標準識別渦核對比Table 3 Comparison of vortex core identifing results based on negative λ2 and positive Q criteria

在平面流中,由負 λ2與 ?u的復特征值和正Q定義的區域是等效的,如考慮平面流的一般速度梯度式

?u的特征 值σ 的特征 方程為σ2+Q=0,其 中Q=?(a2+bc),故σ=±(?Q)0.5,Q>0 時?u為復特征值,由式

負 λ2要求a2+bc<0,即Q>0,故對于平面流這三個定義等價.

2 結果與分析

2.1 基于局部最小壓力標準

此標準通過Burgers 渦流[48]和Lamb-Oseen 渦流[49]兩個簡單例子進行說明,見圖5.軸對稱 Burgers旋渦代表了完整的 N ?S 方程的少數已知精確解之一,渦流由純渦疊加在無旋轉的背景應變流上組成,該流動是不可壓縮的,在柱坐標系中,長度是相對于Burgers 長度尺度的無量綱化,時間是相對于a?1的無量綱化,其中a是背景流場的應變率,其解可以寫為

其中,Γ 為環量,ν 為流體黏度.Lamb-Oseen 渦流的方位速度可寫為

此解對應于 N ?S 方程的自相似解,半徑隨時間增加,有a2(t)=a2(0)+4νt.這里背景流擴散被忽略了,渦流被認為是凍結的,a恒定,黏度的影響由雷諾數來表征.進行無量綱化,取渦旋半徑a=1 和環量 Γ=2π ,角速度和軸向渦度分別可表示為

對于Burgers 渦流,由于沒有旋轉應變率,渦承受軸對稱張力,而張力通過渦的擴散來平衡[50];渦核的長度尺度是固定的,壓力峰值的長度尺度,即?p/?r=0點的半徑(圖5(a)),隨著軸上渦度的增加而增加.Lamb-Oseen 渦流是從初始線渦開始衰減的軸對稱渦流[3](圖5(b)),當半徑r=4(νt)0.5時,相對渦度ω′幾乎消失,而相對壓力p′仍有一個顯著值(最大負值約為0.25).可以看出,渦核和相關局部低壓區域之間存在固有的比例差異,這使得使用局部最小壓力標準進行渦識別成為一個問題.

圖5 相對渦度和壓力分布(實線為渦度,虛線為壓力,黑虛線為ω0=1,向下橙紫綠黃藍 ω0 分別為 2,3,···,6)Fig.5 Relative vorticity and pressure distribution (solid line is vorticity;dotted line is pressure;black dotted line is ω0=1 ;downward orangepurple-green-yellow-blue lines are ω0=2,3,···,6 respectively)

2.2 流線或跡線標準

圖6 中顯示出了在三個不同移動參考系中的軸對稱Lamb-Oseen 渦流的流線型式.使用流線或跡線標準識別渦,將模糊在不同參考系下的渦流,并不可避免地在包含不同速度的多渦平流的湍流中會失敗.

圖6 Lamb-Oseen 渦在不同參考系下的流線(快速點比任何一點的速度都快)Fig.6 Streamline of vortex under different reference systems (the speed at the fast point is faster than that at any point)

2.3 渦度值

螺旋槳伴流中的渦度等值剖面如圖7 所示,可見槳軸向渦度變化很大,具有強大的核心動力學特性[47].盡管實際槳轂區域只有一個連續的渦流柱,但考慮到 | ω| 值的設置具有一定的隨意性,等值面可能沿著渦旋軸終止并被分成多段,有必要設置不明確的足夠低水平 | ω| 值來識別渦,這是渦度值標準識別渦的缺點.

圖7 渦度大小等值面Fig.7 Contours of vorticity magnitude

2.4 張量 T 的復特征值

為了對此標準渦識別方法進行物理解釋,考慮運動參考系中Lamb-Oseen 渦的流線,見圖6(b).在以相同速度移動的粒子A和B的參考系中,渦流核心內的A周圍出現閉合的流線型態(用粗線標記),表示正 ?,而在渦流核心外的B處出現了鞍形,? 為負.

2.5 張量 T 的第二不變量 Q 或運動渦度數 Nk

螺旋槳伴流中的Q=0 剖面如圖8 所示,可見葉尖渦、葉片后緣渦和轂渦都清楚地被識別,也可看出他們之間的相互誘導作用,與圖7 有明顯差別.在螺旋槳下游大約1D后,葉尖渦和后緣渦之間的相互作用變得模糊;在2D之后,轂渦、葉尖渦和后緣渦相遇,進一步模糊尾跡;尾跡近場主要形成一些小半徑渦核,遠場中渦核直徑增大,軸中心的渦核半徑在下游增加得更多.

圖8 渦核識別 Q=0Fig.8 Vortex core identification Q=0

2.6 張量 S 2+?2 的特征值 λ2

螺旋槳伴流中的 λ2=0 剖面如圖9 所示,可以看出其與圖8 的高度一致性,說明采用Q標準和 λ2標準對于螺旋槳伴流中渦識別基本等效.

圖9 渦核識別 λ2=0Fig.9 Vortex core identification λ2=0

3 結論

(1)局部低壓標準比較直觀,但深究考慮其黏性和非定常的影響后,明顯不足;跡線或流線顯然不能滿足伽利略不變性,會使辨別變得混亂;渦度大小需要規定其閾值,具有一定不確定性,且也會識別不是渦的渦片;張量T的復特征值也會有識別不出的區域;速度梯度張量的第二不變量Q標準和對稱張量的第二特征值λ2標準能更好地識別渦核,這兩種標準有時等效.

(2)對于船用螺旋槳的伴流采用Q標準和 λ2標準等效,都能很好地捕獲螺旋槳伴流中的渦核,且能很好地解釋螺旋槳尾流渦脫落的機理.