梢渦空化CFD數值方法

韓寶玉，熊鷹，劉志華

(海軍工程大學 船舶與動力學院，湖北 武漢 430033)

梢渦空化是螺旋槳及高速艇水翼空化的主要組成部分.梢渦空化發生時，不僅會輻射大量噪聲，誘發船體振動，還會使槳后舵產生嚴重剝蝕.在設計中要控制甚或消除梢渦空化的這些不利影響，需要準確了解梢渦流場的物理特征及空化的形成機理.

實驗研究［1-4］在揭示螺旋槳及水翼梢渦空化的形成和發展方面做了大量工作.與實驗相比，數值方法成本低且能得到流場內部多數細節問題，有助于人們理解內部流場的特征，因此應用數值方法模擬升力體梢渦空化非常有意義.

當前，應用RANS方法結合多相流理論模擬空化流形態及其分布，以此預報空化對流體機械性能的影響是數值研究的熱點.RANS方法在計算梢渦流場時的難點是數值離散精度，渦核處網格尺度及湍流模型;而多相流理論中難點在于如何建立空化模型.

Dacles-Mariani等［5］研究了RANS方法中數值離散和網格尺度對梢渦計算精度的影響，認為對方程非粘性項離散要達到五階精度，渦核橫剖面每個方向至少有15個網格節點才能準確計算渦的流場分布.Spall等［6］通過對網格的特殊處理，認為在梢渦流場計算中只要渦核處網格尺度合適，二階離散的精度就可以達到工程要求.在湍流模型研究方面，Churchfield等［7］發現現有的未經過旋轉和曲率修正的一方程和兩方程湍流模型都過高預報了渦核處的湍流粘性耗散，導致渦強迅速減小.而經過旋轉和曲率修正［8］的湍流模型有效抑制了渦核處的粘性耗散，提高了計算精度.

空化模型方面，黃勝等［9］應用完全空化模型模擬了二維翼穩態空泡.韓寶玉等［10］用改進 VOF (volume fraction)模型［11］模擬了NACA66二維翼型的非穩態空化分布，對低空泡數時非穩態空泡的發展、脫落現象進行了研究.

本文以翼剖面為NACA16020翼型的橢圓形水翼為研究對象，采用經過旋轉和曲率修正的代數雷諾應力湍流模型預報了水翼梢渦流場分布，并與實驗結果進行了對比.在此基礎上，用改進VOF模型模擬了不同空泡數下梢渦空化的軌跡和形態，并對梢渦空泡消失空泡數進行了預報.

1 數學方法

1.1 數值方法及湍流模型

控制方程為

梢渦流場中的高曲率及強旋轉流動導致湍動粘度的各向異性影響較大，應用渦粘性湍流模型模擬的流場分布與實驗結果明顯不符.為了克服這些缺點，Wallin等［12］提出了顯式代數雷諾應力模型(EARSM)，該模型對雷諾應力的輸運方程進行了簡化處理，忽略了方程中一些關乎次要流動現象的高階項，建立了雷諾應力與應變率張量和渦張量間的非線性關系式，考慮了湍流各向異性的影響.實驗表明梢渦流場的高旋度及流線的高曲率造成的近似剛體旋轉運動穩定了梢渦內流場分布，使梢渦尾流中的湍動能迅速衰減.為了捕捉渦的這種穩定效應，Wallin等［13］提出了適用于EARSM湍流模型的旋轉和曲率修正方法，以減小渦核內的湍動粘性耗散.韓寶玉等［14］通過對常用兩方程模型和代EARSM模型進行旋轉和曲率修正，研究了湍流模型對機翼梢渦流場數值模擬的影響.發現經過修正的EARSM湍流模型計算結果較其他湍流模型與試驗結果吻合更好.以此為依據，本文計算中選用經過旋轉和曲率修正的EARSM湍流模型.

1.2 空化模型

根據Bakir等［11］提出的理論，混合相由液相、汽相和不可壓離散氣體3部分組成.假定混合相速度相同，混合相密度:

把氣體和液相看作一體(用下標m表示)，假定氣體的質量分數為常數，其密度ρm可表示為

式中:ρ為密度，y為質量分數，下標d表示離散氣體，l表示水，v表示汽相，y*=ρ*/α*/ρ.該混合相的連續方程為

式中:α為體積分數，汽相體積分數αv=1－αm;=－，為混合相和汽相之間質量交換率.顯然只需知道的表達式，聯立式(1)～(5)即可求得汽相體積分數.

分析球形氣泡脈動過程可得

式中:Rb為氣泡半徑，N為單位體積內氣泡的數目.根據Rayleigh-Plesset方程，忽略方程二階項及粘性項，可得

式中:pv、p、ρl分別為氣化壓力、流場壓力、水的密度.而

根據式(6)～(8)可得

式中:Fc是經驗系數，用于修正氣泡膨脹和收縮的速度.在式(4)～(9)中，yd、Rb、Fc均為經驗常數.根據文獻［10］，當pv－p＞0時，Fc=100;當pv－p≤0時，Fc=0.005.計算中認為yd=1.3×10－7，Rb=1×10－6m.

2 計算模型及網格劃分

本文的計算模型是翼剖面為NACA16020翼形的半橢圓形水翼，翼最大弦長C=0.475 m，展長L=0.712 5 m.實驗［1］在空泡水洞中完成，工作段的橫剖面尺寸為1.14 m×1.14 m.來流速度 V= 10 m/s，基于弦長和來流速度的雷諾數為4.75× 106，翼攻角α=10°.

為了驗證計算方法的準確度，本文建立的計算域與實驗一致.原點在翼根弦長中心處，入口和出口與翼導邊和隨邊距離分別為2.5 C，10.0 C.計算域及坐標系如圖1所示.

圖1 計算域Fig.1 Domain used in computations

根據計算域的幾何形狀，本文主要用H型網格對計算域進行整體網格劃分，見圖2.考慮到梢渦的形成與翼邊界層內流動有關，為了較準確地模擬邊界層內流場，本文應用O型網格對翼壁面附近的網格進行了處理，以提高網格質量，第一層網格尺度y+在1～30之間.

圖2 計算域網格示意Fig.2 Overall view of computational domain grid

由于梢渦渦核內徑向速度梯度較大，為了減小該區域的數值離散誤差，本文對梢渦渦核區域的網格進行加密，使渦核內徑向網格節點數達到50×50以上，整個計算域網格總數為500萬.

在邊界條件設置上，來流為了保證速度和雷諾數與實驗一致，設置進口的總壓為151 170 Pa，對應來流速度為10 m/s，入口處的湍流強度為0.15%，渦粘比為5;出口邊界條件為質量流率12 957 kg/s;翼表面及計算域上下表面為不可滑移物面邊界條件.計算在商用CFD軟件CFX中完成.

3 計算結果與分析

本文首先計算了水翼梢渦流場速度分布，通過與實驗對比驗證該方法的可行性.在此基礎上應用相變空化模型對梢渦空化特性進行了研究.

3.1 尾流中梢渦流場特性

尾流中梢渦的近似剛體旋轉運動給梢渦流場的數值預報增加了難度，在模擬梢渦內流場特性時，應用標準代數雷諾應力模型(EARSM)模型和經旋轉和曲率修正后的湍流模型(EARSM-CC)計算結果有明顯差別.

圖3和圖4分別對比了x方向3個位置處梢渦內無因次軸向速度(U/U∞)和切向速度(Vt/U∞)沿z軸分布的計算結果和實驗結果.

圖3 梢渦內軸向速度分布Fig.3 Profiles of velocity in tip vortex

從圖3中可以看出，在各個x位置處，未經修正的代數雷諾應力模型預報的軸向速度峰值明顯低于實驗結果，尤其在x/C=0.3時其峰值低于來流速度，與實驗結果截然相反.這是因為隨梢渦向下游發展，梢渦的近似剛體旋轉加強，使得渦粘性耗散減小，而未經修正的湍流模型過高地預報了梢渦內湍流粘性耗散，導致梢渦內速度遞減過快.而經過旋轉和曲率修正的代數雷諾應力模型則有效減緩了這種湍動粘性耗散，尤其在近場x/C=0.1處與實驗結果吻合很好.而該位置處正是梢渦空化初生的關鍵位置，為后續空化初生的研究提供了保證.

從圖4可以看出，經過旋轉和曲率修正的湍流模型準確預報出了各個x位置處的切向速度分布，與實驗結果吻合較好;而未經修正的湍流模型則過低地預報了切向速度的峰值，與實驗結果差距明顯.從實驗結果可以看出，在尾流近場各個x位置處，當z＜0時，梢渦內流場切向速度存在明顯的波動，且曲線呈現嚴重的不對稱性.文獻［15］認為造成這種波動的原因是水翼后二次渦對梢渦的干擾，計算準確模擬出了二次渦的這種干擾作用.

圖4 梢渦內切向速度分布Fig.4 Profiles of tangential velocity in tip vortex

圖5為應用修正后的湍流模型計算的梢渦渦強無因次軸向分量(ωx×C/U∞)在x/C=0.1和0.2處云圖(視角為從下游往上游看)，圖中可見梢渦主渦和二次渦.從圖中可以看出，在x/C=0.1處二次渦的強度比較大，對主渦的影響較明顯，這從圖4計算結果中可以看到.隨著梢渦往下游發展，二次渦與壁面發生分離，并逐漸融合到主渦中去.這種現象可以從圖6中看到，圖中主渦所指曲面為ωx×C/U∞=100等值面，二次渦所指曲面為 ωx× C/U∞=－20等值面.

圖5 梢渦的形成和卷曲Fig.5 Formation and rollup of the tip vortex

圖6 無因次軸向渦強等值面圖Fig.6 Iso-surfaces of normalized axial vorticity

3.2 梢渦空泡數值模擬

前文計算結果顯示，在計算水翼梢渦近場速度分布時，本文所用的湍流模型及網格尺度是可信的.在此基礎上，本文引入改進VOF模型，對水翼梢渦空泡消失的臨界空泡數及梢渦空泡形態進行了預報研究.

梢渦空化的初生是流體中氣核膨脹直至形成肉眼可見的空泡的過程，由于氣核的膨脹需要一定時間，而相變模型基于空化在瞬間形成和潰滅的假定，因此該模型在預報梢渦空化從無到有的初生過程存在先天缺陷.在實驗測量中空泡的初生空泡數受實驗條件的影響較大，而消失空泡數則比較穩定，因此人們通常用消失空泡數作為空泡發生的衡準數.和空泡初生相比，空泡的潰滅是短暫的，因此應用相變模型預報梢渦空泡的消失空泡數具有一定的可行性.

圖7為水翼在空泡數σ=2.5、2.6、2.7、2.8時汽相體積分數αv=0.1的等值面計算結果.從圖中可以看出，在該空泡數范圍內，水翼只出現梢渦空化，而片狀空泡尚未出現.隨空泡數增大梢渦空泡的橫截面尺寸減小，空泡條紋越來越短.文獻［16］認為當空泡尺寸達到1～3 mm時認為空化消失，本文以2 mm為判斷依據.當空泡數達到2.6～2.7時，該等值面橫截面尺寸達到2 mm左右，此時可定義為空泡消失，其空泡數即為空化消失空泡數σd，計算值與實驗［1］測量結果σd=2.65吻合較好.

圖7 不同空泡數時梢渦水汽含量等值面圖(αv=0.1)Fig.7 Iso-surfaces of water vapour volume fraction αv=0.1 in different cavitation number

4 結論

本文應用RANS方法計算了翼剖面為NACA16020翼型的橢圓形水翼梢渦流場分布，并結合改進VOF相變空化模型對梢渦空泡特性進行了研究，得出以下結論:

1)在計算近似剛性體旋轉運動的梢渦流場時，未經旋轉和曲率修正的湍流模型過高地預報了尾流場中的粘性耗散，導致渦核內軸向速度減小過快.而經過修正的代數雷諾應力湍流模型較好地改進了計算結果，與實驗吻合較好.

2)本文應用的數值方法準確預報出了梢渦內主渦和二次渦的關系，與實驗測量結果相符，為后續對梢渦空泡的研究奠定基礎.

3)應用相變空化模型預報梢渦空泡的軌跡和消失空泡數是可行的，其預報結果與實驗吻合較好.

［1］FRUMAN D H，DUGUE C，CERRUTTI P.Tip vortex rollup and cavitation［C］//19th Symposium on Naval Hydrodynamic.Washington D C，USA，1992:633-654.

［2］ARNDT R E A，ARAKERI V H，HIGUCHI H.Some observations of tip-vortex cavitation［J］.Journal of Fluid Mechanics，1991，229:269-289.

［3］ARNDT R E A，MAINES B H.Viscous effects in tip vortex cavitation and nucleation［C］//20th Symposium on Naval Hydrodynamic.Washington D C，USA，1994:268-289.

［4］MAINES B H，ARNDT R E A.Tip vortex formation and cavitation［J］.Journal of Fluid Engineering，1997，119: 413-419.

［5］DACLES-MARIANI J，KWAK D，ZILLIAC G.On numerical errors and turbulence modeling in tip vortex flow prediction［J］.International Journal for Numerical Methods in Fluids，1999，30:65-82.

［6］SPALL R E.Numerical study of a wing-tip vortex using the Euler equations［J］.Journal of Aircraft，2001，38(1):22-27.

［7］CHURCHFIELD M J，BLAISDELL G A.Near field wingtip vortex computation using the WIND code［R］.Washington DC:AIAA，2006.

［8］SPALART P R，SHUR M.On the sensitization of turbulence models to rotation and curvature［J］.Aerospace Science and Technology，1997，5:297-302.

［9］HUANG Sheng，HE Miao，WANG Chao，et al.Simulation of cavitating flow around a 2-D hydrofoil［J］.Journal of Marine Science and Application，2010，1:63-68.

［10］韓寶玉，熊鷹，陳雙橋.對二維翼空化流動的數值模擬［J］.水動力學研究與進展，2009，24(6):740-746.

HAN Baoyu，XIONG Ying，CHEN Shuangqiao.Numerical simulation of cavitation around 2-dimentional hydrofoil［J］.Chinese Journal of Hydrodynamics，2009，24(6):740-746.

［11］BAKIR F，REY R，GERBER A G.Numerical and experimental investigations of the cavitating behavior of an inducer［J］.International Journal of Rotating Machinery，2004，10:15-25.

［12］WALLIN S，IOHANSSON A.A complete explicit algebraic Reynolds stress model for incompressible and compressible flows［J］.Journal of Fluid Mechanics，2000，403:89-132.

［13］WALLIN S，IOHANSSON A.Modeling streamline curvature effects in explicit algebraic Reynolds stress turbulence models［J］.International Journal of Heat and Fluid Flow，2002，23(5):721-730.

［14］韓寶玉，熊鷹，葉金銘.湍流模型對三維翼梢渦流場數值模擬的影響［J］.航空學報，2010，31(12):2342-2347.

HAN Baoyu，XIONG Ying，YE Jinming.Effect of turbulence models on numerical simulation of tip vortex flow field behind 3D wing［J］.Acta Aeronautica et Astronautica Sinica，2010，31(12):2342-2347.

［15］ARNDT R E A，KELLER A P.Water quality effects on cavitation inception in a trailing vortex［J］.Journal of Fluids Engineering，1992，114:430-438.

［16］盛振邦，劉應中.船舶原理(下)［M］.上海:上海交通大學出版社，2004:68.