胡常莉,王國玉,陳廣豪,王復峰,趙靜
(北京理工大學 機械與車輛學院,北京100081)
當液體內部的局部壓強降低到液體的汽化壓強以下時,在液體內部或液固交界面上就會產生包含蒸汽或氣體的空穴(空泡),這種現象稱為空化。空化現象是航行體水下航行過程最重要的相關流動現象之一[1],它直接影響航行體的水動力特性及操縱穩定性。研究回轉體的空化對工程實際應用具有重要的意義。
多年來,國內外學者對回轉體的空化流動進行了廣泛研究。Rouse[2]等通過實驗研究得到了回轉體表面壓力分布情況,并為數值計算模型的評價提供了可靠的數據。May[3]研究了自然及通氣空化狀態下的回轉體的動力特性并分析了動力脈動及特征頻率與空泡形態之間的關系。劉樺等[4]采用高速攝影技術對1/4 平頭軸對稱體的空泡形態進行了實驗研究,發現空泡形態的斷裂會產生低頻脈動現象。數值計算方面,文獻[5 -6]的研究表明湍流模型對于空化流動的預測有十分重要的影響。Johansen等[7]提出基于標準k-ε 模型和大渦模擬的濾波器模型(FBM),文獻[8]應用FBM 模型計算了繞圓盤空化器的超空化流場,研究了流場結構及水動力特性。Girimaji 等[9]提出了一種混合RANS/LES 的局部時均化湍流模型(PANS),Frendi 等[10]用PANS 和DES模型預測了后臺階跌坎的湍流流動,發現在相同的網格尺寸和邊界條件下PANS 模型的計算結果更精確。文獻[11 -13]構建了不同空化模型源相來描述空化流動中的汽液相間傳輸速率,這種基于質量傳輸方程的空化模型廣泛應用于空化流動數值計算中。
為了深入地了解繞平頭回轉體的空化流動特性,本文綜合運用實驗及數值計算的方法,研究了不同空化數下繞平頭回轉體空穴的發展及空泡的脫落特性,對比分析了不同空化數下空穴的發展過程及脫落細節。
實驗在一循環式空化水洞中進行,如圖1所示,空化水洞系統主要由蓄水池、穩流除氣罐、電機及調速系統、軸流泵、真空發生裝置、實驗段及管路組成。實驗段截面為0.19 m × 0.07 m 的矩形,長度為0.7 m. 通過其上下部及前側面的透明有機玻璃窗觀察空泡形態。實驗中空化數的定義為

式中:p∞、u∞、ρ 和pv分別為距實驗段上游入口210 mm處參考斷面上的平均靜壓強、斷面平均速度、水的密度和汽化壓強。實驗時,保持流速為8.5 m/s,其對應的雷諾數Re=1.7 ×105,通過真空泵調節參考斷面的壓強進而調節空化數。

圖1 空化水洞示意圖Fig.1 Schematic diagram of cavitating tunnel
圖2為高速全流場顯示系統布局示意圖。流動顯示實驗時,采用3 臺1.2 kW 鏑燈照明,高速攝像機記錄流場中空化形態的演變歷程,采集速度設置為5 000 幀/s.

圖2 高速全流場顯示系統布置圖Fig.2 Layout of high-speed video camera system

采用均質平衡流模型,汽液兩相混合物的連續性方程和動量方程:式中:下標i 和j 分別代表坐標方向;ρm、u 和p 分別為混合密度、速度和壓強;μ 和μt分別為混合介質的層流和湍流黏度;ρl、ρv分別為水的密度和水蒸汽的密度。
本文中采用PANS 模型[9]對湍流流動進行求解。PANS 模型的湍動能ku和耗散率εu的輸運方程分別為

式中:Pu、εu分別為湍動能的產生項和耗散項。湍動粘度:

與標準RANS k-ε 兩方程比較,在PANS 模型中,主要對耗散系數C*ε2做出了如下修正:

在高雷諾數的流動中,fε值通常取1. 當fk=1時,說明湍流控制方程復原到RANS 模型;當fk=0時,表示數值計算過程沒有湍流模型的引入,為直接求解的方式。本文通過對現有軟件的二次開發,嵌入了此模型。
空化流動計算中,選用Kubota 空化模型[13],其表達式如下:

式中:RB為簡化氣泡半徑;pv為汽化壓強;ρl、ρv分別為液體密度和蒸汽密度;αnuc為氣核的體積分數;Fe和Fc分別是蒸發和凝結常數項。Kubota 空化模型重點考慮了空化初生和發展時空泡體積變化的影響,適于模擬空化的非定常特性。
計算采用與實驗幾何尺寸相同的模型。計算區域及邊界條件設置如圖3所示,邊界條件設置為速度進口和壓力出口,流動區域的上下左右邊界采用自由滑移壁面條件,平頭回轉體表面采用絕熱、無滑移固壁條件。

圖3 計算域及邊界條件Fig.3 Computational domain and boundary condition
圖4給出了平頭回轉體及其一個縱向切面的網格圖。為了更準確計算空化流動,在平頭回轉體周圍近壁區域進行了網格加密,近壁面y+值為20 ~100,滿足壁面函數要求。回轉體的前端采用O 型結構的網格,這樣可以較好地匹配平頭回轉體頭部的形狀。數值計算的空化數及雷諾數大小均與實驗值保持一致。

圖4 平頭回轉體及其縱切面網格Fig.4 Grid of blunt body and its longitudinal section
圖5(a)和圖5(b)分別給出了σ=0.9 時,由高速錄像和數值計算得到的不同時刻的空穴形態變化。圖5(b)中同時給出了相應的回轉體表面壓力分布云圖。從圖5(a)中可以看出,空穴為由環繞在回轉體的頭部附近小尺度的空泡團組成。在該空化條件下,空穴的發生在位于距肩部約0.2D ~1.5D的范圍內(D 為回轉體直徑);在周向上空穴并沒有完全包裹回轉體。結合圖5(b)可知,聚集狀的空泡團發生在回轉體頭部的低壓區域,且該區域的壓力分布呈現明顯的非定常性。實驗及數值結果均發現,該工況下空穴的發展過程并沒有明顯的周期性,且空泡團的脫落現象表現為尾部小尺度空泡團的瞬時潰滅消失。

圖5 瞬時空穴形態圖(σ=0.9)Fig.5 Evolution process of cavity with time (σ=0.9)
為了說明空化進一步發展階段空穴隨時間的變化情況,圖6和圖7分別給出了σ =0.7 與σ =0.6時,空穴形態在其一個發展周期內的演變過程。相比于σ=0.9 時的空穴形態,發現這兩個工況下的空穴向回轉體的頭部發展并呈附著狀,另外,由于空化數的降低,空化得以進一步發展,空穴形態均以大尺度空泡團的形式包裹著整個回轉體的頭部區域,同時進行著規律、周期性的演變:開始時,空穴呈橢球狀,較穩定地包裹著回轉體頭部;隨后,空穴表面開始凹陷;當t0+ 4.4 ms 時,在距回轉體肩部約0.3D 的位置處,空穴斷裂;然后,空穴尾部出現空泡團逐漸脫落的現象;最后,空穴又覆蓋了回轉體的頭部區域,此時空穴形態便完成了一個周期的變化過程。雖然這兩個工況下的演變過程基本相同,然而仔細對比圖6和圖7則可發現二者存在明顯的差異:其一,當空穴斷裂之后,σ=0.7 時的空穴幾乎完全被一分為二,而σ =0.6 時的空穴在一定時間后又重新融合在了一起;其二,這兩個工況下空穴發展的準周期分別為15.4 ms 和37.4 ms,即σ =0.6 時的發展周期更長些且空穴可以在較長的時間內穩定地包裹著回轉體頭部而不出現大尺度空泡團的脫落現象。

圖6 空穴形態隨時間的演變(σ=0.7)Fig.6 Evolution process of cavity with time (σ=0.7)
空穴的斷裂與空泡的脫落特性是研究非定常空化問題的關鍵。從3.1 節的研究中發現,在不同空化數下,繞平頭回轉體的空泡脫落過程是不同的。為了進一步研究繞平頭回轉體空泡的脫落特性,本節將先后討論σ=0.7 和σ=0.6 時,空泡的脫落過程及其特點。

圖7 空穴形態隨時間的演變(σ=0.6)Fig.7 Evolution process of cavity with time (σ=0.6)
表1中分別給出了σ =0. 7 時,由實驗及數值計算得到的空泡脫落細節,其中,實驗結果共列出了A、B、C 3 組不同的空泡脫落過程,數值計算結果同時給出了回轉體的表面壓力分布云圖。從表1中可以看出,實驗結果及數值計算結果比較一致地展示出了空泡脫落的全過程:空穴斷裂—裂痕前面附著空穴逐漸增長—裂痕后面大尺度空泡團翹起并逐漸脫落。對比由實驗得到的A、B、C 3 組結果可知,空穴斷裂的位置基本相同,即位于距回轉體肩部約為0. 3D 的位置處,且空穴斷裂為U 型渦狀。

表1 空泡脫落細節(σ=0.7)Tab.1 Details of cavity shedding (σ=0.7)
圖8統計了表1中A、B、C 3 組空穴斷裂位置隨時間的變化情況及脫落空泡的運動速度。統計時,始終以斷裂處后面的某一定點為統計對象。3 條位置曲線可以反映出脫落的空泡隨主流的運動情況,結合其速度曲線可知,空泡脫落的速度隨時間是波動變化的,其波動范圍為0 ~1.8u0(u0為來流速度),它們的平均速度約為0.6u0. 由表1中空穴的裂痕隨著時間的推移而逐漸變寬的趨勢可推測出,裂痕前面的附著空穴增長的平均速度小于0.6u0.
類似地,圖9給出了σ =0.6 時,空泡的脫落細節圖。從圖9中可以發現,空穴斷裂的位置也位于距肩部約0.3D 處。當t0+0.8 ms 時,空穴裂痕增大且與回轉體軸線呈一定的角度。隨著時間推移,由于裂痕上游的空穴發展較快,裂痕逐漸模糊,裂痕前后的空泡團逐漸融合在一起。當t0+2.0 ms 時,空穴的下游處又發生空穴斷裂,之后,斷裂后的兩部分空穴再次融合。正是這種“斷裂”又“融合”的作用,抑制了大尺度空泡團的卷起及脫落現象,取而代之的是小尺度空泡團在空穴尾部脫落潰滅。圖10給出了另一組空穴的斷裂及融合的演變過程,從圖10 中可以清楚地看到,與σ =0.7 時相類似,在反向射流的作用下空穴的斷裂亦呈U 型渦狀。

圖8 空泡脫落過程及脫落速度隨時間的變化情況(σ=0.7)Fig.8 Change of cavity shedding and its velocity with time

圖9 空泡脫落細節(σ=0.6)Fig.9 Details of cavity shedding (σ=0.6)

圖10 空穴斷裂與融合(σ=0.6)Fig.10 Cavity shedding and its amalgamation (σ=0.6)
空泡的脫落特性與流場中的反向射流是密切相關的,為了研究反向射流對空穴的作用,圖11給出了空化數σ=0.6 時,在反向射流的作用下空穴的發展情況,同時給出了計算獲得的回轉體的表面壓力分布云圖及其一個縱切面上的速度矢量分布圖。結合空穴形態及回轉體的表面壓力分布情況可知,在空穴閉合區域的后面存在較大的逆壓區域,尤其是最長空穴的后面對應著局部高壓區,此處會首先誘導一股反向射流,隨著時間推移,空穴尾部的高壓區逐漸沿回轉體的周向方向擴展,同時結合速度矢量分布圖可知,反向射流也逐漸增強,導致空穴不斷回縮。由空穴形態圖可以看出,回轉體周向各點的反向射流向前推進的速度是不同的,導致空穴斷裂呈U 型渦狀,這一特點與實驗觀測的結果相一致。

圖11 反向射流推進過程(σ=0.6,上圖為空泡形態與壓力云圖,下圖為速度矢量與壓力云圖)Fig.11 Development of re-entrant jet (σ=0.6,the images above are cavity shapes and pressure contours,the images below are velocity vectors and pressure contours)
本文采用高速錄像方法和數值計算方法研究了繞平頭回轉體的空穴發展及空泡脫落特性,所得結論如下:
1)繞平頭回轉體空穴的發展過程和空泡的脫落細節均受空化數的影響。空化數較大時,空穴的發展及脫落過程沒有準周期性;隨著空化數的減小,空穴的發展及脫落過程具有明顯的周期性。研究發現,空化數越小,空穴的發展周期越長,脫落的空泡尺度越小。
2)繞平頭回轉體的空泡準周期性脫落過程主要包括空穴斷裂與空泡脫落兩個過程。統計發現,空穴斷裂的位置大約位于距回轉體肩部0.3D 處。實驗及數值結果均表明,由于反向射流的作用使空穴的斷裂呈U 型渦狀。
3)在一定的空化條件下,繞平頭回轉體的空化由于斷裂和融合兩個過程的存在,將使空穴呈現更加穩定性的特征。
References)
[1]權曉波,李巖,魏海鵬,等.大攻角下軸對稱航行體空化流動特性試驗研究[J]. 水動力學研究與進展,2008,28(6):662 -667.QUAN Xiao-bo,LI Yan,WEI Hai-peng,et al. An experiment study on cavitation of underwater vehicle’S surface at large angles of attack[J]. Chinese Journal of Hydrodynamics,2008,28(6):662 -667.(in Chinese)
[2]Rouse H,McNown J S. Cavitation and pressure distribution,head forms at zero angel of yaw,studies in engineering[R]. Bulletin:State University of Iowa,1948:32.
[3]May A. Water Entry and the Cavity-Running Behavior of Missiles[R]. Silver Spring:Naval Sea Systems Command Hydroballistics Adcisory Committee,1975.
[4]劉樺,何友聲,趙崗. 軸對稱空泡流的脈動性態研究[J]. 上海力學,1997,18(2):99 -104.LIU Hua,HE You-sheng,ZHAO Gang. A study on pulsation of axisymmetric cavitation flows[J]. Shanghai Journal of Mechanics,1997,18(2):99 -104. (in Chinese)
[5]Wu J,Wang G Y,Shyy W. Time-dependent turbulent cavitating flow computations with interfacial transport and filter based models[J]. International Journal for Numerical Methods in Fluids,2005,49(7):739 -761.
[6]Coutier-Delgosha O,Reboud J L,Delannoy Y. Numerical simulation of the unsteady behavior of cavitating flows[J]. International Journal for Numerical Methods in Fluids,2003,42(5):527-528.
[7]Johansen S T,Wu J,Shyy W. Filter-based unsteady RANS computations [J]. International Journal of Heat and Fluid Flow,2004,25(1):10 -21.
[8]余志毅,王國玉,顧玲燕,等.圓盤空化器超空化繞流流場結構及動力特性的數值分析[J].兵工學報,2008,29(12):1444 -1449.YU Zhi-yi,WANG Guo-yu,GU Ling-yan,et al. Numerical analysis of structure and dynamic characteristics of suppercavitating flow around a disc cavitator[J]. Acta Armamentaii,2008,29(12):1444 -1449. (in Chinese)
[9]Girimaji S,Abdol-Hamid K S. Partially averaged navier stokes model for turbulence:implementation and validation[C]∥The 43rd AIAA Aerospace Sciences Meeting and Exhibit. Reno,Nevada:AIAA,2005:502.
[10]Frendi A,Tosh A,Girimaji S S. Flow past a backward-facing step:Comparison of PANS,DES and URANS results with experiments[J]. International Journal for Computational Methods in Engineering Science and Mechanics,2006,8(1):23 -38.
[11]Kunz R,Boger D,Stinebring D,et al. A preconditioned implicit method for two-phase flows with application to cavitation prediction[J].Comput & Fluids,2000,29(8):849 -875.
[12]Singhal A K,Athavale M M,Li H,et al. Mathematical basis and validation of the full cavitation model[J]. Journal of Fluids Engineering,2002,124(3):617 -624.
[13]Kubota A,Kato H,Yamaguchi H. A new modeling of cavitating flows:a numerical study of unsteady cavitation on a hydrofoil section[J]. Journal of Fluid Mechanics,1992,240(1):59 -96.