基于渦量矩理論的繞振蕩水翼渦動力學分析1)

郝會云 * 劉韻晴 * 魏海鵬 張孟杰 ** 黃 彪 *,

* (北京理工大學機械與車輛學院,北京 100081)

? (北京宇航系統工程研究所,北京 100076)

** (中國北方車輛研究所,北京 100072)

引言

動邊界繞流問題廣泛存在于船舶推進[1-3]、水力發電[4]、海洋工程[5-6]等和國家安全以及國民經濟相關的重要工程領域.振蕩水翼作為軸流泵、水輪機等旋轉水力機械中的重要基礎研究單元,繞其復雜流場結構對水力機械的動力特性有重要影響[7-9].

研究發現繞動態邊界的流場常常存在流態轉捩、流動分離、失速旋渦演化等復雜流動行為.早在20 世紀90 年代,文獻[10]便通過實驗測量了一振蕩NACA0012 機翼表面的壓力分布與渦通量,成功證實了機翼吸力面上動態失速渦的出現、發展與分離過程.根據旋渦出現的不同位置分別定義有前緣渦(LEV)、尾緣渦(TEV)等.以沉浮翼型為研究對象,Ellington 等[11]發現昆蟲和鳥類扇動翅膀時瞬間的高升力歸因于機翼前緣渦的增長.Kissing 等[12]對做俯仰和沉浮組合運動的平板前緣渦的形成和分離行為進行了深入研究,發現到達一定攻角時,二次渦會切斷前緣渦從剪切層積累環量的路徑,隨后前緣渦向下對流并與平板分離.文獻[13]采用歐拉法與拉格朗日擬序結構方法深入分析了繞二維俯仰振蕩翼型的流場演化特征,發現前緣渦誘導較低吸力面上回流從而導致尾緣渦生成.

隨著大量研究深入,人們發現流場演化、旋渦發展與物體的水動力特性存在顯著的關聯.流動分離和動態渦脫落常常導致泵/水輪機在失速狀態下運行,對應時刻流場結構及水動力特性極其復雜,常伴有強烈的振動[14-15].文獻[16]采用實驗和數值方法研究了繞振蕩Clark-Y 水翼的空化流動結構和水動力特性,發現逆時針尾緣渦的形成會導致升力下降,而二次渦(SV)的形成導致水動力載荷的波動.Alam 和Muhammad[17]通過理論分析與數值計算,闡明了振蕩水翼周圍流動結構的演變與產生推/阻力的關系.Xia 和Mohseni[18]通過對比二維振蕩平板的旋渦發展與氣動升力,發現前緣渦的緩慢對流產生的負升力較小,而尾緣渦的快速脫落產生的正升力較大.這兩種升力貢獻之間的差異導致了一個持續約兩個弦長行程的正升力.Tang 等[19]實驗測量了繞NACA66 的水動力特性,發現水翼的低頻振蕩是由前緣分離渦沿流向覆蓋水翼吸力面的長短切換引起的.

在上述研究的基礎上,為了深入了解渦動力學對水動力特性的影響,亟需建立旋渦結構與水動力的定量表征關系.相關研究者對此探索提出了許多著名的理論,如K-J 定理[20-22]、普朗特升力線理論[23]等,然而在實際應用中,這些理論均有一定局限性[24],只能用于穩態流場,無法捕捉流動分離、旋渦、激波等現象.對于在不可壓縮流體中,作任意運動或變形的物體,文獻[25]提出了“渦量矩理論”,揭示了其所受的力和力矩只取決于渦量場的一階及二階矩的變化率,該理論能夠適用于非定常流,但必須求解物體運動產生的整個渦量場(包括啟動渦),除簡單近似外,不能用于有渦量逸出的有限區域.在“渦量矩理論”的啟發下,又陸續發展了其他用渦量表述的流體動力學理論,例如邊界渦量流理論,只需通過物面上動力學量的分布即可求解物體水動力,得到廣泛應用[26-28],其中孫茂[29]基于邊界渦量流理論解釋了昆蟲飛行的物理機制.然而邊界渦量流理論尚不能清晰地揭示流體內部影響空氣/水動力的流動結構.其中文獻[24]應用導數矩變化手段發展得到有限域渦量矩理論,獲得了在任意有限域流體中自由運動物體所受合力的表達式.該公式通過分別選取物面和無窮遠處作為有限域外邊界,能夠復現上述邊界渦量流理論與渦量矩理論,同時解決了渦量矩理論只能適用于無限域的弊端.文獻[30]后續將該理論推廣至高雷諾數非定常流計算,并證實了理論的普適性.

基于以上調研,盡管針對振蕩水翼旋渦演化與水動力特性的研究很多,但是目前尚未有針對水翼局部旋渦結構與其動力特性的綜合量化分析.因此,本工作采用數值模擬方法以及有限域渦量矩理論,結合實驗數據,研究了俯仰振蕩NACA66 水翼的渦動力學特性及其動力響應規律,以期為更好地控制水力機械的動力特性或高效利用水翼運動提供參考數據.

1 研究方法

1.1 控制方程與數值方法

在基于N-S 方程的計算框架內,采用標準k-ωSST 湍流模型[15].其基本控制方程如下:

質量守恒方程

動量守恒方程

式中,下標i和j分別代表坐標方向,u和p分別為速度和壓力,μt為湍流黏性系數,ρ為密度.

由Menter 提出的k-ωSST 湍流模型如下

其中,Pk和Pω為湍流生成項,Dk為湍流耗散項,F1和F2為混合函數,S為剪應力張量的常數項.

計算模型的計算域尺寸、水翼的相對位置以及邊界條件設置均與Ducoin[31]在法國海洋工程研究中心(Research Institute of French Naval Academy)的空化水洞中所完成的實驗保持一致,以實現后續對比.如圖1 所示,邊界分別設置為速度入口和壓力出口,入口速度5 m/s,水翼表面采用絕熱、無滑移的固壁條件.計算域主要包括靜態外域和動態內域兩部分,其中矩形靜態域的長度為16c,高度為1.28c,圓形動態域的直徑為1.1c.靜態域和動態域之間的數據交互由CFX 表達式語言(CEL)控制,并通過設置動態域的振蕩運動來模擬NACA66 水翼的振蕩運動.振蕩速度為6°/s,振蕩幅度為15°,單個周期內振蕩軌跡為0°→15°→0°,定義水翼逆時針向上旋轉過程攻角為α+,水翼順時針向下旋轉過程攻角為α-.用以對比的實驗數據為測得的對應靜態固定攻角下水翼的升力系數.

圖1 計算域與邊界條件設置Fig.1 Computational domain and boundary conditions

圖2 給出了計算域網格劃分情況,在計算網格中使用了結構化網格策略,對水翼近壁區域網格,尤其是水翼前緣(LE)、尾緣(TE)和尾跡區域進行加密,在邊界層設置了500 個節點,使其滿足y+=yuτ/ν≈ 1.

圖2 水翼周圍網格劃分示意圖Fig.2 Computational grids around the hydrofoil

1.2 基于有限域渦量矩理論的合力計算方法

在任意有限域流體中自由運動物體所受合力的具體表達如下[30]

上述表達式中,Vf是以任意控制面 Σ 為外邊界和物面?B為內邊界的流體區域,各瞬時變量均為RANS 模型下直接得到的均化結果,每一項都有明確的物理意義.

式(5)右側第1 項為由翼型運動的非定常慣性效應引起的局部流體渦量變化;第2 項為Lamb矢量或渦力項;第3 項為Lamb 矢量在有限域外邊界處的線積分,能夠捕捉有限域外脫落旋渦結構對水動力的影響;第4 和第5 項分別為物體運動和變形對水動力的作用以及有限域外邊界上的黏性作用.

其中Lamb 矢量積分項代表一個非線性水動力,需要速度u和渦量ω共存,而力的方向總和二者垂直,所以不做功.在局部意義下,這種力正是升力的特點;FB取決于物體自身的運動,而與流動無關,可視為流場的驅動機制;Ft是時均計算模型下的湍流附加項,不影響總合力計算,可忽略該項的求解.對于二維運動模型,上述表達式中n=2,k=1,ρ和μ=ρν分別為流體的密度和黏度,此時物體所受阻力D即為公式各矢量積分項的X方向分量之和Fx,所受升力L為公式各矢量積分項的Y方向分量之和Fy.以Lamb 矢量為例,其阻力方向和升力方向分量分別記作lx和ly.

其中,U∞為無窮遠來流速度,c為水翼弦長,s為水翼翼展長度.

2 結果與討論

2.1 基于Lamb 積分項升力計算分析

流場中的復雜旋渦可由渦量場與流線捕捉.對于二維流場只存在垂直于流動平面方向的渦量,其定義為

為分析繞振蕩NACA66 水翼非定常旋渦結構對其水動力產生的影響,采用基于有限域渦量矩理論的合力計算公式建立局部旋渦與升力的聯系.有限域設置為能夠包括邊界層區域的矩形域,同時為了觀察翼面上方的旋渦演化行為,矩形下游邊界需要超過水翼尾緣.圖3 中 Σ1,Σ2,Σ3為三個任意矩形有限域外邊界,各域右邊界距離水翼尾緣距離分別為ΔX1=0.5c,ΔX2=1.5c,ΔX3=4.17c.圖4 對式(5)在不同域下的計算能力進行了準確性驗證,已知CFX流體計算軟件中采取的內置升阻力計算原理為經典的物面積分

圖3 繞振蕩NACA66 水翼有限域示意圖Fig.3 Diagram of finite domains around the pitching NACA66 hydrofoil

圖4 不同有限域升、阻力系數的計算值Fig.4 Predicted lift/drag coefficient in different finite domains

通過與CFD 輸出結果對比,有效驗證了式(5)在應用中的有效性,同時說明該公式對升阻力的計算與有限域的選取無關.這是由于有限域外的脫落尾渦對物體受力的影響主要體現在流動的非定常性上,如誘導水翼附近渦層的波動等.

圖5 展示了式(5)各分項對升力的影響.圖5(a)給出了數值計算直接得到的升力系數隨水翼轉角的變化趨勢,以及采用式(5)計算的結果,并與同樣來流速度、雷諾數、靜態固定攻角下Ducoin[31]在實驗中獲取的NACA66 水翼的升力系數進行了對比.可見,采用該公式所得到的Cl曲線在幅值和相位上與CFD 計算結果、實驗測量數據均吻合較好,為下文基于算例流場預測的準確性開展定量分析奠定基礎.在圖5(b)中,Lamb 矢量積分項對升力的貢獻最大,幾乎起全部作用,且與CFD 直接結果基本吻合;渦沖量變化項和邊界渦量切割項與渦的脫落、對流行為有關,在動態失速階段由于旋渦不斷生成并向下游對流穿過有限域外邊界,兩項存在較小幅值的波動,其余時間貢獻近似為零;在圖5(c)中,由于低速振蕩且頻率較低,剛體運動項和外邊界黏性作用項提供的升力系數量級極小,對總升力的貢獻可以忽略.

圖5 升力公式 (5)的分解Fig.5 Decompositon of lift Eq.(5)

綜上可認為Lamb 矢量積分項在升力中占據絕對主導地位,下文將采用這一項來探究局部流動結構對升力的貢獻情況.

圖6 對比了不同積分域下Lamb 矢量積分項的計算值與CFD 結果,隨著積分域的增加,Lamb 矢量積分項計算得到的升力系數曲線較CFD 曲線向右偏離更多,即輕微延遲.為保證能夠充分全面地觀察到水翼表面上的旋渦發展,同時兼顧計算的準確性,選擇右側邊界距離水翼尾緣0.5 倍弦長的矩形域作為積分域,即上述 Σ1所圍繞的域,開展下面的工作.

圖6 不同有限域下Lamb 矢量積分項的升力計算能力Fig.6 Predictive ability of the Lamb vector integral in different finite domains

2.2 旋渦動力學特性定量分析

為深入分析旋渦結構與動力特征之間的相互作用,結合流場速度流線與渦量分布云圖以及上述升力演化曲線,在圖7 中將水翼振蕩的整個周期分為三個階段:線性增加段stage I (α+=0°～ 14°)、失速波動段stage II (α+=14°～ 12°)和線性減小段stage III(α-=12°～ 0°).

圖7 流場階段劃分Fig.7 Stage division

在線性增加段有以下幾個特征階段.

(1)α+=0°～ 4°:在該階段中,水翼的升力系數隨著攻角的增大呈線性增長的趨勢.此時繞振蕩水翼流場呈現準層流狀態,如圖8(a)所示.

(2)α+=4°～ 6°:升力系數在該階段有一明顯的拐點.這是因為繞水翼前緣流場出現了層流向湍流的過渡轉捩(laminar separation bubble,LSB).如圖8(b)所示,當α2=4.08°時,水翼前緣處無層流分離現象產生,在水翼尾緣X/C=0.7 位置處,由于逆壓梯度的作用,形成LSB 區域.如圖8(c),當水翼旋轉到α3+=5.4°時,層流向湍流的過渡轉捩區域逐漸轉移至水翼前緣位置,這一過程導致升力曲線出現拐點.

圖8 線性增加段典型時刻流場Fig.8 Typical flow fields in linear increasing stage

(3)α+=6°～ 14°:升力系數在此階段不斷增大.在水翼吸力面上,出現一順時針的尾緣渦-TEV,如圖8(d)所示;在水翼向上旋轉的過程中,-TEV 逐漸增大,并向水翼前緣發展;α5+=13.4°時,如圖8(e)所示,在水翼尾緣的吸力面上,-TEV 的下游出現一逆時針的尾緣渦 +TEV.

圖9 給出了水翼在向上振蕩過程中,分離點及再附著位置隨轉角的變化趨勢,其位置根據Prandtl流動分離判據給定 (?ux/?y)y=0=0.從圖中可以看出,繞振蕩水翼流場呈現準層流狀態時,層流分離點位于水翼尾緣X/C=0.6～ 0.7 之間,并隨攻角增加緩慢前移;當水翼攻角位于α+=4°～ 6°時,層流分離點由水翼尾緣迅速移動至前緣,并在靠近前緣頂點處發生層流向湍流的過渡轉捩;此后全局流場轉變為湍流流動,特別是當進入動態失速階段時,流場中的大尺度分離渦不斷生成、相互作用并對流脫落,而分離點始終位于接近水翼頂點位置.

圖9 分離點及再附著位置隨攻角的變化情況Fig.9 Evolution of separation points and reattachment positions with the angle of attack

在失速波動段(α+=14°～α-=12°),升力系數曲線存在三次升力峰值,對應著三次周期性的旋渦發展與脫落,每個周期內存在相似的旋渦演化行為,波動峰值分別位于α+=13.92°,α-=14.73°,α-=13.51°,現以旋渦發展的第二個周期為例,如圖10 所示,研究流場結構演化與水動力特征.

圖10 失速波動段典型時刻流場Fig.10 Typical flow fields in dynamic stall stage

(1)α+=14.76°～α-=14.73°:該階段水翼升力系數處于上升階段,水翼吸力面在-TEV 和 +TEV 的基礎上誘導形成一個順時針前緣渦-LEV,如圖10(a)所示,此時繞振蕩水翼同時存在三個旋渦,隨后,-LEV 和-TEV 同時向水翼中央移動,直至攻角α10=14.73°時融合形成一個新的-LEV,幾乎覆蓋水翼整個吸力面,如圖10(b)所示,此時升力系數達到局部峰值.

(2)α-=14.73°～ 14.32°:升力系數隨攻角減小呈近似線性下降趨勢.這個過程中,-LEV 向水翼尾緣快速移動,并與 +TEV 相互作用直至二者完全脫落.如圖10(c)所示,此時-LEV 即將完全脫離翼面,升力達到局部最小值.

由此推斷-LEV 的充分發展以及-LEV 與+TEV的相互作用脫落是導致升力系數曲線波動的重要影響因素.

線性減小段與線性增加段的流場變化規律相似.

(1)α-=12°～ 5°:該階段升力系數隨攻角減小而線性減小,最初在水翼尾緣同時存在-TEV 和 +TEV,如圖11(a),隨后-TEV 和 +TEV 相互作用,+TEV 脫落,如圖11(b)所示,水翼吸力面僅剩下-TEV.

圖11 線性減小段典型時刻流場Fig.11 Typical flow fields in linear decreasing stage

(2)α-=5°～ 4 °:該階段流場行為近似于線性增加段,繞振蕩水翼的層流分離點由水翼前緣向尾緣移動,導致升力曲線出現拐點.

(3)α-=4°～ 0°:升力系數隨攻角減小線性下降,如圖11(c),繞振蕩水翼流場由湍流再度轉變為層流.

由上述分析可知,繞振蕩NACA66 水翼流場中存在著復雜的非定常旋渦,尤其是在動態失速階段,旋渦演化呈現出與升力的強關聯性,從流場觀測上來看,振蕩水翼獲得高升力的機制是前緣渦不斷發展并附著在水翼吸力面上,翼后緣的渦層不斷地脫瀉到尾流中去,但是其中的定量關系尚不明確.下面通過對不同旋渦內ly積分來計算其對瞬時升力的貢獻情況.

如圖12 所示黃色封閉虛線為典型時刻下,基于渦量與流線所識別的各個旋渦結構積分域.

圖12 典型時刻局部旋渦結構積分域示意圖Fig.12 Diagram of vortices’ domain of integration at a typical moment

圖13 展示了動態失速階段九個特征時刻的旋渦貢獻情況.可以看到,在升力抵達峰值之前的上升段,存在典型的三渦流場,結合圖10(a),此時-LEV與-TEV 同時向水翼中央移動,+TEV 附著在水翼尾緣.這個過程中-LEV 和-TEV 提供了部分正升力.在α6,α9,α12三個時刻中,-LEV 對總升力的貢獻分別為17%,9%,17%,-TEV 對總升力的貢獻為29%,10%,12%,而+TEV 提供了近似為零的極小負升力;在α7,α10,α11三個時刻,瞬時升力達到局部峰值,如圖10(b),-LEV 幾乎覆蓋整個水翼吸力面,其對瞬時升力貢獻分別為50%,37%,37%;在α12,α13,α14時刻,瞬時升力達到局部最低值,此時-LEV 和+TEV 經過相互作用即將或者完全脫離水翼表面,-LEV 的貢獻接近于0,+TEV 提供了負升力,分別占瞬時升力的-20%,-22%,-44%.

圖13 失速波動段關鍵旋渦貢獻Fig.13 Key vortices’ contributions during the dynamic stall stage

綜上所述,附著的-LEV 和-TEV 提供正升力,+TEV 提供負升力;峰值時-LEV 的貢獻最大接近50%,當升力達到局部最小值時,-LEV 的貢獻約為0,而 +TEV 的負貢獻可高達40%.

上文分析了有限域內不同類型旋渦對升力的凈貢獻情況,然而值得注意的是,旋渦內部不同局部區域所提供的升力有正有負.如圖14 所示為ly云圖,左側為全局流場圖,右側為各個旋渦區域的局部放大圖.對比圖10 發現ly不僅能夠精確捕捉到對應的旋渦區域,同時能夠直觀表達出不同流場結構的升力貢獻,紅色表示對應的旋渦結構提供了正升力,藍色表示對應的旋渦結構提供了負升力.且各個旋渦內部均被分割成多個部分,這是由于流場速度的方向改變導致,而具體旋渦的凈貢獻是區域積分的結果.

圖14 旋渦結構與Lamb 矢量Y 向分量云圖Fig.14 Vortices and ly contours

圖15 追蹤了第二次峰值所對應的旋渦發展脫落周期中,特定-LEV,-TEV,+TEV 對瞬時升力的貢獻大小與演化情況.可以看到,各旋渦所提供的升力與其對瞬時升力的貢獻率保持同步,當融合后的-LEVly積分值達到最大值時,它對瞬時升力的貢獻也最大,接近50%;在發展過程中,+TEV 始終提供負升力,-LEV 與-TEV 在融合之前提供正升力,融合后新-LEV 繼續向下對流,與 +TEV 相互作用脫離水翼表面的過程中,在α-=14.56°時產生短暫的負升力,而這期間+TEV 的負貢獻較大,在α-=14.56°時高于40%.因此失速行為發生時可認為是+TEV 通過誘導有利-LEV 的脫落及自身貢獻負升力兩種途徑從而使升力快速下降.

圖15 一個失速周期內的旋渦升力貢獻演化Fig.15 Contributions of vortices during one stall period

為了進一步探究上述旋渦逸出有限域后將如何影響水翼的動力特性,圖16 給出了在失速階段典型時刻尾流渦分布,包括渦量場和ly云圖.可以看到順時針尾流渦各區域均提供正升力,而逆時針渦各區域均提供負升力.

圖16 尾流旋渦結構與Lamb 矢量Y 向分量云圖Fig.16 Wake vortices and ly contours

結合圖5 中邊界渦量切割項可準確捕捉域外旋渦的貢獻,在α-=14.38°時,域外尾渦提供了負升力,隨后脫落的順時針旋渦(提供正升力)與脫落的逆時針旋渦(提供負升力)不斷逸出有限域而進入域外尾渦區域,由于不斷逸出的逆時針旋渦“攜帶”更高的正升力,大于同一時刻進入域外的順時針渦所提供的負升力,因此域外尾渦對升力的貢獻逐漸上升.在α-=13.96°時,尾渦整體提供了正升力.但總體來看,域外尾渦對升力的貢獻較小,僅在失速階段,域外尾渦提供的升力發生小幅波動,體現了流動的非定常性.

3 結論

(1) 繞振蕩NACA66 水翼流場存在邊界層轉捩、流動分離和動態渦失速等行為,隨著水翼攻角的增加,繞水翼的層流向湍流轉捩點由尾緣移至前緣發生,隨后在水翼尾緣吸力面上出現順時針尾緣渦-TEV,-TEV 不斷增大并向前緣發展,與水翼前緣生成的順時針旋渦-LEV 融合形成覆蓋水翼整個吸力面的-LEV,此后與尾緣吸力面逆時針旋渦+TEV 相互作用直至完全脫落分離.在水翼向下旋轉階段,水翼近壁面流動逐漸由湍流轉變為層流.水翼振蕩過程中,前緣渦與尾緣渦的相互作用引起前緣渦完全脫落分離,直接導致了動力失速現象的發生.

(2)動態失速渦與升力演化密切相關.基于有限域渦量矩理論的升力公式由渦量變化項、Lamb 矢量積分項、邊界渦量切割項、剛體運動項以及外邊界黏性作用項5 項構成,其中Lamb 矢量積分項占主導,能夠用來表達振蕩水翼升力演化,進而捕捉局部旋渦結構對升力的貢獻,并直觀展示不同流動結構對升力的作用方向及大小,結果顯示附著的-LEV和-TEV 提供正升力,當-LEV 發展覆蓋整個吸力面時,其對升力的局部貢獻最大且接近50%,+TEV 提供負升力.然而有限域內任一旋渦內部的不同位置提供的升力有正有負;域外脫落尾渦卻只提供一種貢獻,順時針渦提供正升力,逆時針渦提供負升力.