基于THINC/QQ 格式的超空泡魚雷三維流場特性及水動力性能優化研究

肖 騰,謝彬,2,杜艷平

(1.上海交通大學 海洋工程國家重點實驗室,上海,200240;2.上海交通大學 船海工程數值試驗中心,上海,200240;3.上海交通大學 中英國際低碳學院,上海,201306)

0 引言

自然超空泡是高速航行器周圍的流體因壓力降低至飽和蒸汽壓以下發生相變后形成的較大蒸汽與氣體的空腔。當空泡尺度擴張到足以包括整個航行器時,就形成超空泡。借助超空泡技術,蘇聯研制出了“暴風”(Shkval)超空泡魚雷,使航速提高至200 kn 以上,受到國內外學者的廣泛關注。資料表明[1],水下航行器在有幾乎完全包覆表面的穩定超空泡的條件下,僅有頭部空化器直接與水體接觸,航行器所受到的粘性阻力與壓差阻力顯著降低,相較于全濕流動,總阻力可降低1～2 個數量級。

對于超空泡的發展與機理,許多學者都做了相關的研究工作。Ahn 等[2]進行了多組水洞實驗,系統地分析了超空泡的形態與幾種常用空化器頭型的關系。Saranjam[3]對水下移動物體的超空泡流動進行了數值仿真和實驗研究,數值預測與實驗結果高度一致。賈立平等[4]進行了自然超空泡高速彈射實驗,研究了超空泡的形成過程和空化器參數對超空泡臨界空化數的影響,并分析了超空泡形態與空化器的關系。王瑞等[5]使用高速水洞實驗平臺研究了4 種不同形狀空化器對繞回轉體超空泡特性的影響,給出了相同工況下空化器頭型與空泡尺寸和阻力特性的關系。

然而,超空泡實驗裝置耗資較高,且實驗中的介入式裝置會對空泡流場產生嚴重干擾。隨著空化理論和計算方法的不斷發展,數值仿真方法逐漸成為超空泡問題的重要研究手段。齊江輝等[6]基于雷諾平均類(Reynolds-Averaged Navier-Stokes,RANS)模型對無尾翼超空泡魚雷模型開展了超空泡流仿真,分析了不同空化數和各種圓盤空化器的改型對超空泡形態和減阻效果的影響。Kim 等[7]通過改變空化器形狀和開角,仿真研究了直航超空泡魚雷的空化器優化設計方法,給出一定工況下的空化器最優解并進行了拖曳實驗驗證。李雨田[8]通過仿真計算,得到了空泡外形隨空化數和攻角的變化趨勢,分析了航行器流體動力特性變化,并給出了影響其變化的相關因素。趙瀟雨等[9]使用FLUENT 軟件對4 種典型空化器頭型進行了超空泡數值仿真研究,分析了自然超空泡尺寸與空化器形狀、尺寸和空化數之間的關系。除了改變空化器形狀,在材料表面涂覆親疏水性涂層可實現流動減阻。黃超等[10]實驗研究了表面親疏水性對小球入水后空泡的影響,其結果顯示,超疏水表面的小球由于表面液膜很快產生流動分離,相比于親水表面的小球更易于空泡的產生。不過,表面親疏水性對魚雷超空泡的影響還有待進一步研究。

不同于此前基于商業軟件的數值仿真研究,文中基于THINC/QQ 方法,利用剪切應力傳輸(shear stress transport,SST)k-ω湍流模型和Schnerr-Sauer空化模型在開源軟件OPENFOAM 上開發了基于非結構網格的不可壓縮空化流求解器,對超空泡魚雷的三維流場特征進行了數值仿真研究,并用經驗公式對仿真結果的可靠性進行了驗證。在此基礎上分析討論了不同空化數、空化器形狀以及魚雷表面浸潤性對超空泡特征和水動力性能的影響規律。

1 數值模型

1.1 多相流控制方程

文中基于流體體積(volume of fluid,VOF)函數多相流模型,假設將整個空化流動流場中的氣液兩相視為一種混合均質,共享速度場、壓力場和湍流參數,統一求解Navier-Stokes 方程

式中:u為速度;p為修正壓力;μ為流體的動力粘度;φ為液相流體所占的體積分數;為氣液兩相間的質量傳遞率;F=σκ?φ為表面張力項,其中σ=0.072為表面張力系數,κ=?·為單位向量的散度,n=?φ為體積分數的梯度。

汽液兩相通過體積分數值平均得到混合均質的粘度μ和密度ρ,且

式中:ρl和 μl分別為水的密度和粘度;ρv和 μv分別為水蒸氣的密度和粘度。

1.2 VOF 自由界面捕捉方法

為了更好地捕捉超空泡形態,該研究采用VOF高精度算法THINC/QQ[11]。該方法采用雙曲正切函數重構網格單元內的自由界面,對于三維情況,雙曲正切函數為

式中,β是用于控制界面過渡區域厚度的參數。界面方程Pi(x,y,z)+di=0為2 次多項式

式中:cstr(s,t,r=0,1,2,s+t+r≤2)為控制界面形狀的參數,由界面的單位法向量和曲率確定;(xc,yc,zc)為 網格的重心。系數cstr確定后,式(6)中僅有確定界面位置的系數di未知,可根據目標網格的VOF函數均值求解

式中:雙曲正切函數H的網格積分平均值可以用數值積分點近似;系數di通過Newton-Raphson 方法計算,然后求解以下方程更新VOF 函數

式中,下標iup代表邊界面上游方向的網格單元。

1.3 湍流模型

由于超空泡流動具有很高的雷諾數和非定常特性,所以文中引入了SSTk-ω湍流模型求解RANS方程。SSTk-ω模型是由Menter[12]提出,結合了kω模型和k-ε模型的近壁面和遠壁面區域,被廣泛地用于超空泡流動仿真中。該模型中的湍動能k和比動能耗散率ω的計算方法如下

式中:Γ,σk,σω,β,β*是模型常數;混合函數F1為

渦黏系數定義為

式中:a1=0.31;混合函數F2定義為

1.4 空化模型

由于空化流動相比全濕流動存在著相變現象,要對其中的質量傳遞過程進行準確仿真,就要引入空化模型。在OPENFOAM 中,式(2)中的質量傳遞源項被分解為蒸發和液化過程的分量

式中,和m˙-分別為液化和蒸發過程的質量傳遞率。

此處使用了Schnerr-Sauer 空化模型[13],其質量傳遞源項是由Rayleigh-Plesset 方程推導而來,其源項方程如下

式中:Cc和Cv為方程系數,默認推薦值均為1;φnuc為數值成核項,用以在計算中初始化空化過程;psat為飽和蒸汽壓;RB為氣泡半徑。

為了兼顧計算精度和效率,文中采用具有時空2 階精度的有限體積法對控制方程進行離散。其中時間推進采用2 階龍格-庫塔方法,對流項和擴散項分別用Van Leer 格式(迎風格式)和線性格式,速度-壓力耦合則采用PISO(pressure implicit with splitting of operator)算法。時間步長設為自適應,滿足最大庫朗數為0.3 的條件。

2 計算模型

2.1 超空泡魚雷外形

文中使用的計算模型參考Serebryakov 給出的無尾翼魚雷模型[14]。該模型是超空泡航行器數值仿真研究中被廣泛使用的模型之一。文中的超空泡魚雷模型由前部空化器、錐體頭段、平行中段和尾部噴管組成,其尺寸根據“暴風”(Shkval)超空泡魚雷做了優化,旨在給出更貼近真實尺度的數值仿真結果,其總長3.28 m,直徑0.2 m,圓盤空化器直徑0.08 m,具體的外形和尺寸參數如圖1 所示。

圖1 超空泡魚雷模型外形尺寸示意圖Fig.1 Shape and dimensions of the supercavity torpedo model

2.2 計算域與網格劃分

計算域及網格示意圖如圖2 所示。計算域總長度為60D,高度和寬度為20D,近壁區域進行網格加密,向流場邊界逐漸稀疏。邊界條件設置如下:入口為速度入口條件,流速按照空化數(σ=(p-psat)/(0.5ρlu2))來設定;出口為壓力出口條件,按照水深10 m 設置;魚雷表面為壁面條件。研究表明[15],航行器形成自然超空泡的條件是空化數σ ≤0.1,選取45～95 m/s 范圍內的多組航速,對應空化數范圍為0.0965～0.2160 進行計算。如表1 所示,計算前使用多組不同疏密程度的網格進行網格收斂性測試。表1 對比了不同計算網格得到的空化數 σ=0.0270條件下魚雷模型的超空泡尺度數據。對比結果顯示,較為粗糙的網格1 與網格2 得到的數值結果收斂性差,與其他網格的差距較大;網格4 較網格3 進一步加密,得到的超空泡尺度結果變化很小。在確保結果收斂性的基礎上兼顧計算效率,選取網格3 的174 萬網格用以后續計算。

圖2 計算域及網格劃分示意圖Fig.2 Diagram of computational domain and mesh

表1 超空泡尺度網格收斂性測試Table 1 Mesh convergence test on the supercavity scale

2.3 浸潤性邊界條件

超疏水表面特征是通過設置魚雷表面的邊界條件實現的。如圖3 所示,基于超疏水表面的宏觀滑移特點,使用部分滑移邊界條件。根據接觸角給定滑移長度,滑移速度us與滑移長度Ls的關系式為

圖3 部分滑移邊界條件Fig.3 Boundary conditions of partial slip

式中,Ls根據Zhang[16]的統計數據確定,接觸角θ=160°時Ls=2.2。接觸角是衡量材料親疏水性的重要指標(見圖4)。

圖4 接觸角與親疏水性的關系Fig.4 Relation between contact angle and hydrophobicity

2.4 計算模型驗證

文中選取了1 組經典錐形空化器的案例進行模型驗證。在同樣的工況下,文中數值仿真空泡結果與Javadpour 等[17]的水洞實驗結果的對比如圖5 所示。圖中給出了2 組不同速度下的穩定空泡形態,每組上方是文中的數值仿真結果,下方是實驗結果。從圖中可以看出,在不同速度對應的不同空化數條件下,該計算模型對空泡尺度的變化規律預測與實驗結果一致。對空泡長度進行測量和比較發現,數值仿真得到的空泡長度與實驗結果相差5%以內。文中的數值仿真結果與實驗結果的吻合度較高,說明文中的計算模型可以很好地對超空泡流動進行仿真。

圖5 空化器后空泡形態數值仿真結果與實驗結果對比Fig.5 Comparison of the cavity behind the cavitator between experimental results and the numerical results

3 仿真結果分析

3.1 自然超空泡演化過程

圖6～圖8 分別給出了3 組不同空化數工況下的空泡形態演化過程。圖6 為 σ=0.0646工況下空泡形態演化,對應流速為55 m/s,此時魚雷模型并不能產生完整的超空泡,頭部空泡發展至肩部后即為此空化數下的穩定狀態。圖7 為σ=0.0399工況下空泡形態演化,對應流速為70 m/s,此時魚雷模型可以產生長度與模型接近的超空泡,模型的超空泡臨界空化數即σ ≈0.04。圖8 為σ=0.0270工況下典型的自然超空泡演化過程,由圖可見,自然超空泡在演化過程中呈現以下特點:在空化初始階段,僅有空化器末端和尾段末端存在規模較小的空泡;隨著空化數降低,錐體頭端和平行中段連接處的肩部也產生了肩部空泡;3 處主要空泡隨著空化數降低長度將繼續增加,直至互相融合;融合后的超空泡體積將繼續增長,直至穩定。

圖6 σ=0.0646工況下超空泡形態演化過程Fig.6 Evolution process of supercavity shape at σ=0.0646

圖7 σ=0.0399工況下超空泡形態演化過程Fig.7 Evolution process of supercavity shape at σ=0.0399

圖8 σ=0.0270工況下超空泡形態演化過程Fig.8 Evolution process of supercavity shape at σ=0.0270

值得注意的是,在前方空泡與后方空泡融合的過程中,會發生后方空泡體積縮小、部分潰滅的現象。圖9 為空化數 σ=0.0399時頭部空泡與肩部空泡融合時肩部空泡及壓力分布圖。此時,頭部空泡閉合區末端存在高壓區,隨著頭部空泡的發展,此高壓區逐漸靠近肩部空泡,導致其局部液化,體積收縮,直至與前方空泡完全融合后繼續向后發展。主空泡發展至尾段附近時末端的高壓區同樣會在與尾部空泡融合過程中抑制其發展。

圖9 局部空泡融合位置壓力分布Fig.9 Pressure distribution at the position of partial cavity fusion

圖10 給出了空化數 σ=0.0399時臨界超空泡流動過程中魚雷的水動力系數隨時間發展的曲線。從圖10(a)的阻力系數CD曲線中可以看出,在超空泡形成的初期,航行器所受阻力較高,根據數值仿真初始條件及數值結果,空泡在前2 ms 迅速發展,峰值可視為魚雷模型在當前航速下的理想全濕流動阻力;隨著超空泡的發展,航行器所受阻力顯著降低;直至穩定的超空泡流形成后,航行器的阻力系數也降低至一個穩定值,相比超空泡產生初期的峰值減少了88%。圖10(b)中的升力系數CL在超空泡的初生期,由于空化器后部、彈體肩部和尾部有著局部空泡的產生,在此類擾動和近壁區域流場的綜合影響下有著短暫的波動,但幅度較小;在超空泡的發展趨于穩定之后,無尾翼航行器升力值幾乎為零,不再受升力影響。

圖10 超空泡航行器水動力系數時歷曲線Fig.10 Time-history curves of hydrodynamic coefficients of supercavity vehicle

3.2 超空泡形態驗證

為了驗證文中所使用的求解器在超空泡數值仿真方面的準確性,將所得的空泡特征長度和特征直徑(其中超空泡尺度量取方式如圖11 所示)與Reichardt 和Logvinovich 經驗公式[18-19]進行對比驗證。這些經驗公式都是由勢流理論推導而來,適用于小空化數下的超空泡形態計算。

圖11 超空泡尺度測量Fig.11 Measurement of supercavity scale

Reichardt 經驗公式主要將空泡特征尺度與空化數和阻力系數聯系起來,即

式中:lc和dc分別為空泡的長度和直徑,其量取方式如圖11 所示;CD為阻力系數,由以下半經驗公式計算得到

式中,CD0是 σ=0時的阻力系數,根據實驗數據推算得出,根據Self 和Ripken[20]的推導,對于文中的圓盤空化器,取CD0=0.80。

Logvinovich 經驗公式也尤其適用于小空化數,特征尺度與空化數的關系為

式中,a和k為經驗常數,文中分別取值為2.0 和0.95。

取空化數 σ=0.0270為例,將文中的超空泡形態的數值仿真結果與上述經驗公式進行比較,結果如表2 所示,其中相對誤差R和相對誤差L分別對應與Reichardt 公式和Logvinovich 公式的對比結果。如表所示,文中數值仿真結果中的超空泡尺度與經驗公式結果較為接近,驗證了文中數值方法對魚雷模型超空泡尺度預測的可靠性。其中,空泡特征長度相對誤差分別為3.32%和2.56%,空泡特征直徑相對誤差分別為3.17%和4.08%,誤差較低。誤差來源主要是基于勢流理論的經驗公式沒有考慮流體粘性。

表2 超空泡尺度數值結果驗證Table 2 Validation of numerical results of supercavity scale

3.3 空化器形狀尺寸對超空泡的影響

圓盤形是空化器最常見的形狀,其形狀簡單,效果明顯。此外,還給出幾種不同形狀尺寸的空化器改型。各型空化器外形如圖12 所示,其中(a)型為前文所使用的D=0.08 m的圓盤空化器;(b)型為D=0.12 m的圓盤空化器;(c)型和(d)型為錐度不同的D=0.08 m的錐頭空化器;(e)和(f)型為球面半徑分別為 0.04 m 和 0.05 m,截面基圓直徑都是D=0.08 m的球面空化器。

圖12 不同空化器頭型幾何模型Fig.12 Models of cavitators with different head shapes

表3 給出了航速u=90 m/s,空化數σ=0.0241時不同頭型空化器模型在超空泡穩定階段的阻力系數值。由于外形平滑,沒有尖端高壓區,使用小球面半徑球面空化器模型在超空泡階段所受阻力最小,使用小錐角高錐度的錐體空化器次之;使用2 個不同直徑的圓盤空化器模型所受阻力最高,且阻力值大小與空化器直徑成正比。

表3 不同空化器超空泡階段阻力系數Table 3 Drag coefficients for different cavitators at the supercavity stage

使用不同頭型空化器模型的超空泡輪廓如圖13 所示。(b)型較大直徑的圓盤型空化器能夠產生最大體積的超空泡;而(e)型半球面型空化器由于邊緣切向角度不能更好地分離流動,產生的超空泡最細也最短。綜合表3 的阻力系數結果,超空泡的體積尺寸與魚雷模型所受阻力成負相關關系。

圖13 不同空化器產生的超空泡輪廓Fig.13 Contours of supercavity generated by different cavitators

總的來說,扁平比和直徑較大的空化器更容易產生更大體積的穩定超空泡,但是魚雷模型在超空泡穩定階段所受的阻力也更大;使用長徑比和錐度較大、直徑較小的空化器的魚雷模型在超空泡階段所受的阻力更小,但是會在一定程度上損失超空泡的體積。從超空泡減阻的角度考慮,后一類空化器選型方向更利于超空泡魚雷的高速減阻。在實際應用中,需要平衡空化體積與航行阻力之間的關系,空泡體積過小會導致空泡無法完整包覆航行器,從而不能達到減阻效果。

3.4 表面親疏水性對超空泡的影響

為研究表面浸潤性對超空泡的影響,文中選取了上文(e)型半球空化器的魚雷模型,分別設定彈體及空化器表面材質具有超疏水性質(接觸角θ=160°)和超親水性質(接觸角 θ=20°)。

圖14～圖16 分別給出了不同空化數工況下的超疏水表面與親水表面的魚雷模型的阻力系數以及組分里的粘性阻力系數CV曲線。在空化數σ=0.0965工況下,魚雷模型并不能產生完整包覆表面的超空泡,而僅在魚雷頭部產生周期性脈動的局部空泡(見圖6),也造成了圖14(a)中50～60 ms 階段的阻力值波動。超疏水表面整體可降低航行阻力30%左右。而且由圖14(b)可見,全彈體覆蓋的超疏水表面較親水表面在此工況下可以降低粘性阻力約80%。在空化數 σ=0.0399工況下,超空泡能夠完整包覆魚雷表面,總阻力特別是粘性阻力的下降更為明顯。

圖14 σ=0.0965工況下不同接觸角魚雷阻力系數曲線Fig.14 Drag coefficient curves of torpedoes with different contact angles at σ=0.0965

圖15 σ=0.0399工況下不同接觸角魚雷阻力系數曲線Fig.15 Drag coefficient curves of torpedoes with different contact angles at σ=0.0399

圖16 σ=0.0305工況下不同接觸角魚雷阻力系數曲線Fig.16 Drag coefficient curves of torpedoes with different contact angles at σ=0.0305

隨著航速增加,在空化數 σ=0.0305工況下,在空化初期的前20 ms 內,超疏水表面可以降低阻力20%以上。隨著空泡體積擴大,魚雷模型所受阻力中粘性阻力占比減小,阻力差值穩定在10%左右。由圖15(b)的阻力系數曲線可見,超疏水表面從流動初期就保持了很低的粘性阻力水平。在超空泡流動開始階段,超疏水表面模型所受粘性阻力較親水表面模型低約90%;隨著空泡發展,親水表面模型與水接觸面積減小,粘性阻力也逐漸減小。綜合圖14～圖16 來看,隨著航速提高,空化數降低,產生的空泡包覆魚雷彈體的面積也更大,超疏水表面減阻效果會有所下降。從減阻幅值的角度衡量超疏水表面魚雷,控制航速與壓力條件以達到臨界空化數的效果更為明顯。

圖17 為分別使用2 種親疏水表面的(e)型空化器魚雷模型在超空泡初期的超空泡形態。圖17(1)～(5)分別是空泡從空化器末端初生到發展長大過程中5 個時刻對應的形態。超疏水表面魚雷模型產生的空泡長度均略長于同時刻的親水表面魚雷模型,即超疏水表面魚雷模型可以更快地產生空泡,這一點對于實際應用十分有利。根據相關數值結果,超疏水表面魚雷模型還可以產生更大體積的空泡,其頭部空泡長度較親水表面魚雷長約5%,這是因為超疏水的表面特性,超空泡更易附著在彈體表面(見圖17(3)～(5)),而親水表面魚雷模型與空泡間存在一層水膜(見圖17(1)～(2))。隨著空泡發展,由于其閉合區末端的逆壓力梯度產生的回射流將持續存在于彈體與空泡之間,受到表面疏水性的影響,超疏水魚雷表面產生的回射流長度遠小于親水表面魚雷所產生的(見圖17(5)),這有利于超空泡的穩定發展和航行器的持續減阻。圖18 給出了圖17(5)時刻2 種魚雷的速度場。從圖18(b)可見,親水表面魚雷的彈體肩部由于壁面效應存在速度更低的區域,肩部空泡的體積相較于疏水表面魚雷更小。

圖17 2 種魚雷在超空泡初期不同接觸角下的頭部空泡形態Fig.17 Frontal cavity shapes of two kinds of torpedoes with different contact angles at the initial supercavity stage

圖18 2 種魚雷在超空泡初期不同接觸角下的速度場Fig.18 Velocity field of two kinds of torpedoes with different contact angles at the initial supercavity stage

考慮到工程應用實際情況,還研究了超疏水表面材料在魚雷表面不同部位使用的減阻和空泡優化效果的區別。如圖19 所示,將(e)型空化器魚雷模型在軸向劃分為4 個部分:空化器(a)、前段(b)、中段(c)和尾段(d),分別單獨使用超疏水表面進行計算。

一是水資源供需矛盾突出。我國水資源總量位居世界第六，但人均水資源量僅為世界平均水平的1/4。經濟的持續高速或較高增長，人口增長及人均消費水平的持續提高，使得供需矛盾突出，給水安全帶來了日益沉重的壓力，且這種壓力在相當長時間內是不可逆轉的。

圖19 使用超疏水表面的位置(紅色為超疏水表面)Fig.19 Positions of superhydrophobic surfaces employed(red color indicates superhydrophobic surfaces)

圖20～圖21 分別給出了4 種在不同位置布置超疏水表面的魚雷在空化數 σ=0.0305工況下的阻力系數和粘性阻力系數曲線。在彈頭與中段的結合區域布置超疏水表面減阻效果最明顯,尾段次之;在頭部和錐體前端布置超疏水表面減阻效果最不明顯。隨著超空泡的發展,阻力值差距縮小。在中段結合區域布置超疏水表面對不同阻力成分的降低效果最明顯;圖21 中10 ms 后尾段對粘性阻力相較中段降低更明顯,是因為肩部空泡覆蓋了一定面積的彈體。空化器與錐體前端布置超疏水表面對空化初期阻力的降低效果相似,考慮到表面積的差別,在空化器上使用超疏水表面相比錐體前段綜合效果更好。

圖20 不同超疏水位置魚雷阻力系數曲線Fig.20 Drag coefficient curves for torpedoes with superhydrophobic surfaces at different positions

圖21 不同超疏水位置魚雷粘性阻力系數曲線Fig.21 Viscous drag coefficient curves for torpedoes with superhydrophobic surfaces at different positions

4 結束語

基于THINC/QQ 方法,SSTk-ω模型以及Schnerr-Sauer 空化模型,在OPENFOAM 平臺建立三維空化流動求解器,經過與實驗結果對比驗證后,對無尾翼魚雷模型在不同速度、空化器形狀以及彈體表面浸潤性工況下開展了三維超空泡流動數值仿真研究和水動力性能分析,獲得以下結論。

1) 隨著空化數不斷降低,魚雷超空泡的演化過程有著明顯的規律,在空化數達到超空泡臨界值 σ=0.04時,穩定包覆魚雷彈體的超空泡可以使航行阻力下降88%左右。從自然超空泡的演化機理的角度考慮,適當提高航速以保持更穩定和更大體積的超空泡有利于降低航行阻力。

2) 大長徑比、高錐度的空化器有利于超空泡魚雷的減阻,但是會損失超空泡體積。從空化器選型的角度考慮,在減阻以提高航速的同時,也要保證超空泡的尺寸。

3) 超疏水表面可以明顯降低超空泡魚雷在空化初期的航行阻力,尤其是粘性阻力成分,特別是在臨界空泡數工況下。同時,超疏水表面可以加快初始空泡的產生,空泡發展過程中更易于附著在超疏水表面,抑制了空泡內液膜和回射流,有利于提高超空泡的穩定性。

后續的研究工作將基于文中的數值模型,拓展研究通氣超空泡與超疏水技術的應用。