鄔 明,孫善春
(中國船舶重工集團公司第七一○研究所,湖北 宜昌 443003)
超空泡技術是一種可以使水下航行體高效減阻技術,其原理是在水與航行體表面形成穩定的氣層,達到減小水的粘性阻力,從而提高水下航行體的運行速度。水下超空泡航行體所采用的操縱面,要兼顧全沾濕狀態和超空泡狀態時的流體動力布局要求:航行體處于全沾濕狀態時,航行器流體動力布局與常規水下航行體相似,需依靠操縱面提供穩定力矩和操縱力矩。而在超空泡狀態時要盡量減小其所受阻力,同時操縱面必須穿透空泡壁進入水中才能提供所需力矩,該狀態下的操縱面應具有較高的升阻比,因此對于操縱面的研究具有重要意義。
隨著計算性能的快速提高和計算數學理論的不斷發展完善,數值方法已經成為研究問題和解決問題的重要手段之一,它不僅在理論研究領域得到普遍應用,而且在工程實際中也被廣泛使用。求解雷諾平均N-S方程式是當前數值計算的主要方法,需補充湍流模型對方程進行封閉。湍流模型對于數值模擬空化流場的精度具有決定性影響。當前應用最普遍的湍流模型是k-ε 湍流模型[1-2],然而該模型有以下不足之處:在模型中湍流尺度未知;僅限于湍流邊界層壓力相對穩定的情況;其壁面函數在邊界層的修正中難以彌補計算模型與實際現象之間的差距。相對于k-ε湍流模型的以上不足之處,SSTk-ω湍流模型[3]卻擁有以下優點:該模型能夠適應逆壓梯度變化的各種物理現象;可應用于粘性內層,通過對壁面函數的應用,能夠精確模擬邊界層的現象,無需使用較容易失真的粘性衰減函數。本文將SSTk-ω湍流模型應用于操縱面空化流場的仿真分析中,取得了較好的結果。
SSTk-ω湍流模型以下簡稱SST湍流模型。

式中:ρ表示密度;U表示速度;t為時間。

其中

μt和αl分別表示流體動力粘度和水相的體積分數。
若分別以φ1、φ2、φ3表示k-ω 模型、k-ε 模型和SST湍流模型中的函數關系,則SST湍流模型可表示為:

SST湍流模型考慮到湍流剪切應力的輸運,不但能夠對各種來流進行準確的預測,還能在各種壓力梯度下精確模擬分離現象,其綜合了近壁面k-ω模型的穩定性及邊界層外部k-ε模型獨立性的優點,它的計算模擬性能優于后兩者。各系數取值為β′=0.09,α1=5/9,β1=0.075,σk1=2,σω1=2,α2=0.44,β2=0.0828,σk2=1,ωω2=1.168,各個數據的取值取自參考文獻[4]。
選擇Singhal空化模型[5],其中 Rayleigh Plesset方程是控制汽化產生和空泡凝結的理論基礎,該方程描述了液體中空泡的變化過程:

式中r為空泡半徑,Pc為汽化壓強,σvl是汽液兩相間的表面張力系數,質量傳輸率為:

式中,ρg表示飽和蒸汽密度;sgn表示符號函數;NB代表單位體積的汽泡數;F代表經驗因子,汽化過程和凝結過程取值不同,Fvap=50,Fcond=0.01。
針對水下高速航行體在不同運動階段的流體動力特性,對構成操縱面的二維基本翼型進行了相應研究——選用典型的小頂角楔形,經驗表明以上區域邊界對流場的計算結果影響很小。
計算網格的好壞直接影響到數值計算的可行性、收斂性以及計算精度。前處理軟件ICEM CFD是一款成熟的網格劃分軟件,它向用戶提供業界領先的高質量網格技術,其強大的網格劃分功能可以滿足CFD仿真計算的嚴格要求。
對于本文中的問題若采用非結構網格,生成的空泡邊界會出現凹凸不平的毛刺,影響對仿真結果的分析。若想消除毛刺現象,就要對網格進行細分。因此采用基于六面體結構化網格劃分方法,對計算區域進行全六面體結構網格劃分,并利用與湍流強度相關的Yplus對所建模型進行考核。楔形操縱面流場網格圖及局部放大圖如圖1和圖2所示。

圖1 計算區域與網格劃分Fig.1 The computational domain and meshing

圖2 操縱面局部放大圖Fig.2 Zooming in images locally in control surface
邊界條件設定:上游邊界給定來流速度和各流體體積分數;下游邊界給定壓強;流域外邊界為光滑物面邊界。初始條件:設定初始來流速度和各流體體積分數。
采用有限元的有限體積法(FVM),N-S方程中的對流項采用CFX高精度格式,它利用調整混合因子提高計算精度且平均收斂精度。以一階迎風為例,對流項:

式中:φi表示變量φ第i個節點值;φu為變量φ迎風節點值;ξ為混合因子;▽φ為迎風節點的節點梯度;從迎風節點到第i個節點的矢量。
針對不同條件下,ξ會作出相應的調整來滿足計算要求獲得精確的解,例如在第i節點臨近迎風節點的節點梯度變化劇烈,則混合因子趨近于0,反之趨近于1。
根據相關實驗文獻[6],楔形操縱面的典型空泡形態分為三種,如圖3。

圖3 楔形操縱面的典型空化形態(全空化、部分空化和底部空化)Fig.3 Typical characteristics of cavities of control surface(entirely、substantially and bottom enveloped by cavities)
從圖3中不難發現,當來流攻角小于楔形半頂角時,僅在翼型底部產生空化現象,而在上下翼面完全處于全沾濕狀態。此時所產生流體動力和力矩與攻角近似成線性關系,且較為穩定。當來流攻角大于楔形半頂角時,上翼面處于超空泡狀態,而下翼面為全沾濕狀態,此時,同樣可產生穩定的流體動力和力矩,且與攻角成線性關系。當來流攻角與楔形半頂角相當時,除在翼型底部產生空化外,上翼面還將出現間歇性片狀空泡,此時翼型所產生的流體動力及力矩呈非定常脈動分布。
空化數是衡量流體空化程度的物理量,計算公式為:

其中P∞、PV、ρ、V∞分別為某基準位置上的絕對靜壓力、空泡內的絕對靜壓力、流體密度、相對于物體的未擾動流體的速度。研究中設定初始條件P∞、25℃時水的飽和蒸汽壓、V∞和ρ。

L為翼型弦長,升、阻力系數基于速度坐標系。
表1給出了繞楔形操縱面的超空化仿真與試驗結果的對比。

表1 繞楔形操縱面的超空化仿真與試驗結果Table 1 Results of simulation and test of control surface
圖4給出了計算得到的阻力系數和升力系數隨攻角的變化仿真曲線及實驗曲線,兩者基本一致,從中可以看出當保持自然空泡數不變,僅改變來流攻角時,操縱面的阻力系數呈逐漸增大趨勢,而升力系數在攻角增至一定數值后,其增大幅度減緩。隨著攻角的增大,空泡不斷增大,并形成超空泡,使得操縱面全部處于空泡的包裹之中,此時,操縱面的控制效率也將相應地降低。圖5中的升阻比曲線分布表明,當攻角增至一定數值后,升阻比將有所減小,這與臨界攻角和局部空泡的形成有關。

圖4 阻力系數和升力系數隨攻角的變化曲線Fig.4 The change of drag characteristics and lift characteristics by angles of attack

圖5 升阻比隨攻角的變化曲線Fig.5 The change of ratio of lift over drag by angles of attack
CFX后處理可以更加直觀地觀察空泡形態,如圖6所示,不同攻角下,操縱面的密度云圖。從圖中可以看出,當來流攻角很小時(α<3°),楔形操縱面底部產生空泡;當攻角較小時(3°≤α≤3.5°),在操縱面上端面產生局部空泡;而在攻角繼續增大時(α≥3.5°),操縱面完全被空泡包裹,并且呈增大趨勢。這與理論結果可較好地吻合。

圖6 不同攻角下的密度分布圖Fig.6 The contour of the charge density in different angles of attack
圖7 給出的是空化數為0.323時,操縱面的壓力系數分布,可看出,其下翼面上壓力系數分布呈單調遞減趨勢,當攻角a=3°時,上翼面出現局部空泡的壓力系數分布形態,增大來流攻角后,上翼面的壓力系數將降至一恒定值,即上翼面已處于超空化狀態。

圖7 不同攻角,操縱面表面的壓力分布曲線Fig.7 The curve of pressure distribution for control surface in different angles of attack
圖8 給出了不同空化數下的空泡輪廓圖。從圖中可以發現隨著空化數的降低空泡尾部的回射流強度在減弱,這是因為空泡長細比增大,閉合時空泡邊界入射角減小的緣故。這和文獻[7]中的理論分析取得了很好地一致性。

圖8 不同空化數下的空泡輪廓Fig.8 The profile of cavities in different citation numbers
采用SST湍流模型,文章對典型二維操縱面在不同來流空化數及攻角條件下的升阻力特性、空泡形態、表面壓力分布規律進行了詳細的研究,并與相關實驗做了分析比較,結果表明:
(1)減小空化數后,相同攻角條件下的操縱面的阻力系數和升力系數都有所下降。隨著攻角的增大,由于上翼面處于超空化狀態,其升阻比曲線將有所減小;
(2)運用N-S方程及SST湍流模型對均質流場求解,并利用Rayleigh-Plesset空泡模型對二維操縱面的空化流動進行數值模擬,與相關實驗相對比,取得了較好的一致性,證明了所選湍流模型的可行性,為以后的超空泡實驗研究提供參考依據;
(3)SST湍流模型能夠對操縱面的空化流場進行較為準確的模擬,尤其在空泡流場、空泡外形以及尾部回射流方面,可以準確地預測空化流場,從而指導超空泡航行器的外形設計,顯示出了優越性,具有廣闊的應用前景。
[1] 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.
[2] LEROUX J B,COUTIER-DELGOSHA O,ASTO-LFI J A.A joint experiment and numerical study of mechanisms associated to instability of partial cavitationon two-dimensionalhydrofoil[J].Phys.Fluids,2005,17(5):052-101.
[3] COUTIER-DELGOSHA O.Numerical prediction of cavitation flow on a two-dimensional symmetrical hydrofoil and comparison to experiments[J].Journal of Fluids Engineering,2007,129(3):279-291.
[4] KUBOTA A,KATO H.A new modeling of cavitating flows:a numerical study of unsteady cavitation on a hydrofoil section[J].FluidMech.,1992,240(1):59-96.
[5] SINGHAL A K,ATHAVALE M M,et al.Mathematical basis and validation of the full cavitation model[J].ASME Journal of Fluids Engineering,2002,124:617-624.
[6] GARABEDIAN P R.Calculation of axially symmetric cavities and jets[J].Pac.J.Math.,1956,6(4):611-684.
[7] SEMENENKOV N.Artificial supercavitation.physics and calculation[A].RTO AVT lecture series on supercavitating flows[C].Von Karman Institute Brussels,Belgium:2001,206237.
[8] 張宇文.空化理論與應用[M].西安:西北工業大學出版社,2007.