韓寶玉，熊鷹，劉志華

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

梢渦空化是螺旋槳及高速艇水翼空化的主要組成部分.梢渦空化發(fā)生時(shí)，不僅會(huì)輻射大量噪聲，誘發(fā)船體振動(dòng)，還會(huì)使槳后舵產(chǎn)生嚴(yán)重剝蝕.在設(shè)計(jì)中要控制甚或消除梢渦空化的這些不利影響，需要準(zhǔn)確了解梢渦流場(chǎng)的物理特征及空化的形成機(jī)理.

實(shí)驗(yàn)研究［1-4］在揭示螺旋槳及水翼梢渦空化的形成和發(fā)展方面做了大量工作.與實(shí)驗(yàn)相比，數(shù)值方法成本低且能得到流場(chǎng)內(nèi)部多數(shù)細(xì)節(jié)問(wèn)題，有助于人們理解內(nèi)部流場(chǎng)的特征，因此應(yīng)用數(shù)值方法模擬升力體梢渦空化非常有意義.

當(dāng)前，應(yīng)用RANS方法結(jié)合多相流理論模擬空化流形態(tài)及其分布，以此預(yù)報(bào)空化對(duì)流體機(jī)械性能的影響是數(shù)值研究的熱點(diǎn).RANS方法在計(jì)算梢渦流場(chǎng)時(shí)的難點(diǎn)是數(shù)值離散精度，渦核處網(wǎng)格尺度及湍流模型;而多相流理論中難點(diǎn)在于如何建立空化模型.

Dacles-Mariani等［5］研究了RANS方法中數(shù)值離散和網(wǎng)格尺度對(duì)梢渦計(jì)算精度的影響，認(rèn)為對(duì)方程非粘性項(xiàng)離散要達(dá)到五階精度，渦核橫剖面每個(gè)方向至少有15個(gè)網(wǎng)格節(jié)點(diǎn)才能準(zhǔn)確計(jì)算渦的流場(chǎng)分布.Spall等［6］通過(guò)對(duì)網(wǎng)格的特殊處理，認(rèn)為在梢渦流場(chǎng)計(jì)算中只要渦核處網(wǎng)格尺度合適，二階離散的精度就可以達(dá)到工程要求.在湍流模型研究方面，Churchfield等［7］發(fā)現(xiàn)現(xiàn)有的未經(jīng)過(guò)旋轉(zhuǎn)和曲率修正的一方程和兩方程湍流模型都過(guò)高預(yù)報(bào)了渦核處的湍流粘性耗散，導(dǎo)致渦強(qiáng)迅速減小.而經(jīng)過(guò)旋轉(zhuǎn)和曲率修正［8］的湍流模型有效抑制了渦核處的粘性耗散，提高了計(jì)算精度.

空化模型方面，黃勝等［9］應(yīng)用完全空化模型模擬了二維翼穩(wěn)態(tài)空泡.韓寶玉等［10］用改進(jìn) VOF (volume fraction)模型［11］模擬了NACA66二維翼型的非穩(wěn)態(tài)空化分布，對(duì)低空泡數(shù)時(shí)非穩(wěn)態(tài)空泡的發(fā)展、脫落現(xiàn)象進(jìn)行了研究.

本文以翼剖面為NACA16020翼型的橢圓形水翼為研究對(duì)象，采用經(jīng)過(guò)旋轉(zhuǎn)和曲率修正的代數(shù)雷諾應(yīng)力湍流模型預(yù)報(bào)了水翼梢渦流場(chǎng)分布，并與實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比.在此基礎(chǔ)上，用改進(jìn)VOF模型模擬了不同空泡數(shù)下梢渦空化的軌跡和形態(tài)，并對(duì)梢渦空泡消失空泡數(shù)進(jìn)行了預(yù)報(bào).

1 數(shù)學(xué)方法

1.1 數(shù)值方法及湍流模型

控制方程為

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

1.2 空化模型

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

把氣體和液相看作一體(用下標(biāo)m表示)，假定氣體的質(zhì)量分?jǐn)?shù)為常數(shù)，其密度ρm可表示為

式中:ρ為密度，y為質(zhì)量分?jǐn)?shù)，下標(biāo)d表示離散氣體，l表示水，v表示汽相，y*=ρ*/α*/ρ.該混合相的連續(xù)方程為

式中:α為體積分?jǐn)?shù)，汽相體積分?jǐn)?shù)αv=1－αm;=－，為混合相和汽相之間質(zhì)量交換率.顯然只需知道的表達(dá)式，聯(lián)立式(1)～(5)即可求得汽相體積分?jǐn)?shù).

分析球形氣泡脈動(dòng)過(guò)程可得

式中:Rb為氣泡半徑，N為單位體積內(nèi)氣泡的數(shù)目.根據(jù)Rayleigh-Plesset方程，忽略方程二階項(xiàng)及粘性項(xiàng)，可得

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

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

式中:Fc是經(jīng)驗(yàn)系數(shù)，用于修正氣泡膨脹和收縮的速度.在式(4)～(9)中，yd、Rb、Fc均為經(jīng)驗(yàn)常數(shù).根據(jù)文獻(xiàn)［10］，當(dāng)pv－p＞0時(shí)，F(xiàn)c=100;當(dāng)pv－p≤0時(shí)，F(xiàn)c=0.005.計(jì)算中認(rèn)為yd=1.3×10－7，Rb=1×10－6m.

2 計(jì)算模型及網(wǎng)格劃分

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

為了驗(yàn)證計(jì)算方法的準(zhǔn)確度，本文建立的計(jì)算域與實(shí)驗(yàn)一致.原點(diǎn)在翼根弦長(zhǎng)中心處，入口和出口與翼導(dǎo)邊和隨邊距離分別為2.5 C，10.0 C.計(jì)算域及坐標(biāo)系如圖1所示.

圖1 計(jì)算域Fig.1 Domain used in computations

根據(jù)計(jì)算域的幾何形狀，本文主要用H型網(wǎng)格對(duì)計(jì)算域進(jìn)行整體網(wǎng)格劃分，見(jiàn)圖2.考慮到梢渦的形成與翼邊界層內(nèi)流動(dòng)有關(guān)，為了較準(zhǔn)確地模擬邊界層內(nèi)流場(chǎng)，本文應(yīng)用O型網(wǎng)格對(duì)翼壁面附近的網(wǎng)格進(jìn)行了處理，以提高網(wǎng)格質(zhì)量，第一層網(wǎng)格尺度y+在1～30之間.

圖2 計(jì)算域網(wǎng)格示意Fig.2 Overall view of computational domain grid

由于梢渦渦核內(nèi)徑向速度梯度較大，為了減小該區(qū)域的數(shù)值離散誤差，本文對(duì)梢渦渦核區(qū)域的網(wǎng)格進(jìn)行加密，使渦核內(nèi)徑向網(wǎng)格節(jié)點(diǎn)數(shù)達(dá)到50×50以上，整個(gè)計(jì)算域網(wǎng)格總數(shù)為500萬(wàn).

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

3 計(jì)算結(jié)果與分析

本文首先計(jì)算了水翼梢渦流場(chǎng)速度分布，通過(guò)與實(shí)驗(yàn)對(duì)比驗(yàn)證該方法的可行性.在此基礎(chǔ)上應(yīng)用相變空化模型對(duì)梢渦空化特性進(jìn)行了研究.

3.1 尾流中梢渦流場(chǎng)特性

尾流中梢渦的近似剛體旋轉(zhuǎn)運(yùn)動(dòng)給梢渦流場(chǎng)的數(shù)值預(yù)報(bào)增加了難度，在模擬梢渦內(nèi)流場(chǎng)特性時(shí)，應(yīng)用標(biāo)準(zhǔn)代數(shù)雷諾應(yīng)力模型(EARSM)模型和經(jīng)旋轉(zhuǎn)和曲率修正后的湍流模型(EARSM-CC)計(jì)算結(jié)果有明顯差別.

圖3和圖4分別對(duì)比了x方向3個(gè)位置處梢渦內(nèi)無(wú)因次軸向速度(U/U∞)和切向速度(Vt/U∞)沿z軸分布的計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果.

圖3 梢渦內(nèi)軸向速度分布Fig.3 Profiles of velocity in tip vortex

從圖3中可以看出，在各個(gè)x位置處，未經(jīng)修正的代數(shù)雷諾應(yīng)力模型預(yù)報(bào)的軸向速度峰值明顯低于實(shí)驗(yàn)結(jié)果，尤其在x/C=0.3時(shí)其峰值低于來(lái)流速度，與實(shí)驗(yàn)結(jié)果截然相反.這是因?yàn)殡S梢渦向下游發(fā)展，梢渦的近似剛體旋轉(zhuǎn)加強(qiáng)，使得渦粘性耗散減小，而未經(jīng)修正的湍流模型過(guò)高地預(yù)報(bào)了梢渦內(nèi)湍流粘性耗散，導(dǎo)致梢渦內(nèi)速度遞減過(guò)快.而經(jīng)過(guò)旋轉(zhuǎn)和曲率修正的代數(shù)雷諾應(yīng)力模型則有效減緩了這種湍動(dòng)粘性耗散，尤其在近場(chǎng)x/C=0.1處與實(shí)驗(yàn)結(jié)果吻合很好.而該位置處正是梢渦空化初生的關(guān)鍵位置，為后續(xù)空化初生的研究提供了保證.

從圖4可以看出，經(jīng)過(guò)旋轉(zhuǎn)和曲率修正的湍流模型準(zhǔn)確預(yù)報(bào)出了各個(gè)x位置處的切向速度分布，與實(shí)驗(yàn)結(jié)果吻合較好;而未經(jīng)修正的湍流模型則過(guò)低地預(yù)報(bào)了切向速度的峰值，與實(shí)驗(yàn)結(jié)果差距明顯.從實(shí)驗(yàn)結(jié)果可以看出，在尾流近場(chǎng)各個(gè)x位置處，當(dāng)z＜0時(shí)，梢渦內(nèi)流場(chǎng)切向速度存在明顯的波動(dòng)，且曲線呈現(xiàn)嚴(yán)重的不對(duì)稱(chēng)性.文獻(xiàn)［15］認(rèn)為造成這種波動(dòng)的原因是水翼后二次渦對(duì)梢渦的干擾，計(jì)算準(zhǔn)確模擬出了二次渦的這種干擾作用.

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

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

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

圖6 無(wú)因次軸向渦強(qiáng)等值面圖Fig.6 Iso-surfaces of normalized axial vorticity

3.2 梢渦空泡數(shù)值模擬

前文計(jì)算結(jié)果顯示，在計(jì)算水翼梢渦近場(chǎng)速度分布時(shí)，本文所用的湍流模型及網(wǎng)格尺度是可信的.在此基礎(chǔ)上，本文引入改進(jìn)VOF模型，對(duì)水翼梢渦空泡消失的臨界空泡數(shù)及梢渦空泡形態(tài)進(jìn)行了預(yù)報(bào)研究.

梢渦空化的初生是流體中氣核膨脹直至形成肉眼可見(jiàn)的空泡的過(guò)程，由于氣核的膨脹需要一定時(shí)間，而相變模型基于空化在瞬間形成和潰滅的假定，因此該模型在預(yù)報(bào)梢渦空化從無(wú)到有的初生過(guò)程存在先天缺陷.在實(shí)驗(yàn)測(cè)量中空泡的初生空泡數(shù)受實(shí)驗(yàn)條件的影響較大，而消失空泡數(shù)則比較穩(wěn)定，因此人們通常用消失空泡數(shù)作為空泡發(fā)生的衡準(zhǔn)數(shù).和空泡初生相比，空泡的潰滅是短暫的，因此應(yīng)用相變模型預(yù)報(bào)梢渦空泡的消失空泡數(shù)具有一定的可行性.

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

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

4 結(jié)論

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

1)在計(jì)算近似剛性體旋轉(zhuǎn)運(yùn)動(dòng)的梢渦流場(chǎng)時(shí)，未經(jīng)旋轉(zhuǎn)和曲率修正的湍流模型過(guò)高地預(yù)報(bào)了尾流場(chǎng)中的粘性耗散，導(dǎo)致渦核內(nèi)軸向速度減小過(guò)快.而經(jīng)過(guò)修正的代數(shù)雷諾應(yīng)力湍流模型較好地改進(jìn)了計(jì)算結(jié)果，與實(shí)驗(yàn)吻合較好.

2)本文應(yīng)用的數(shù)值方法準(zhǔn)確預(yù)報(bào)出了梢渦內(nèi)主渦和二次渦的關(guān)系，與實(shí)驗(yàn)測(cè)量結(jié)果相符，為后續(xù)對(duì)梢渦空泡的研究奠定基礎(chǔ).

3)應(yīng)用相變空化模型預(yù)報(bào)梢渦空泡的軌跡和消失空泡數(shù)是可行的，其預(yù)報(bào)結(jié)果與實(shí)驗(yàn)吻合較好.

［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］韓寶玉，熊鷹，陳雙橋.對(duì)二維翼空化流動(dòng)的數(shù)值模擬［J］.水動(dòng)力學(xué)研究與進(jìn)展，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］韓寶玉，熊鷹，葉金銘.湍流模型對(duì)三維翼梢渦流場(chǎng)數(shù)值模擬的影響［J］.航空學(xué)報(bào)，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］盛振邦，劉應(yīng)中.船舶原理(下)［M］.上海:上海交通大學(xué)出版社，2004:68.