三維水翼梢渦流場(chǎng)數(shù)值研究

蒲汲君，熊鷹

海軍工程大學(xué)艦船工程系，湖北武漢430033

三維水翼梢渦流場(chǎng)數(shù)值研究

蒲汲君，熊鷹

海軍工程大學(xué)艦船工程系，湖北武漢430033

針對(duì)三維水翼梢渦流場(chǎng)和梢渦空泡初生的問題，分別建立k-ω，DES和LES模型，對(duì)水翼的梢渦流場(chǎng)進(jìn)行計(jì)算研究。為減少誤差，在網(wǎng)格的處理上對(duì)梢渦流域進(jìn)行局部加密，對(duì)未發(fā)生空化時(shí)梢渦內(nèi)的軸向速度和切向速度進(jìn)行計(jì)算，發(fā)現(xiàn)LES模型的計(jì)算結(jié)果與實(shí)驗(yàn)值吻合較好。在此基礎(chǔ)上，討論尾渦的卷曲對(duì)梢渦壓力場(chǎng)的影響，提出了使用氣泡靜力平衡方程計(jì)算初始梢渦空泡數(shù)的方法。

梢渦；氣泡靜力平衡方程；空泡初生；水翼

0 引 言

自1859年Besant［1］的早期研究以來，梢渦空泡一直是學(xué)者們研究的焦點(diǎn)。當(dāng)一個(gè)氣核進(jìn)入到渦核當(dāng)中，如果這里的壓力比臨界壓力更小，那么氣核將會(huì)迅速長(zhǎng)大并形成一個(gè)肉眼可見的空泡?？张菰跐邕^程中會(huì)釋放強(qiáng)烈的輻射噪聲并造成螺旋槳和船舵的剝蝕?？张莸念愋桶ㄉ覝u空泡、片空泡、云空泡和游移型空化。大量實(shí)驗(yàn)發(fā)現(xiàn)，梢渦渦核中心處的局部壓力最低，除了梢部有卸載處理的三維水翼和螺旋槳外，梢渦空泡在所有類型的空泡中最早出現(xiàn)（Keller［2］在其相關(guān)實(shí)驗(yàn)中已有發(fā)現(xiàn)），它產(chǎn)生的噪聲會(huì)對(duì)艦船隱身性產(chǎn)生巨大影響。為抑制梢渦空泡的產(chǎn)生，很有必要準(zhǔn)確計(jì)算出梢渦流域的流場(chǎng)特性。但是，由于梢渦本身的不穩(wěn)定性，準(zhǔn)確把握梢渦流場(chǎng)的特性還是一個(gè)較為困難的問題。解決該問題的有效途徑是進(jìn)行高精度的數(shù)值計(jì)算，它能提供梢渦流域流場(chǎng)的具體細(xì)節(jié)。

Fruman和Arndt等［3-4］對(duì)梢渦的流場(chǎng)分布進(jìn)行了大量研究。與實(shí)驗(yàn)相比，數(shù)值計(jì)算的效率高、成本低，且更容易得到梢渦的流場(chǎng)分布細(xì)節(jié)。因此，使用數(shù)值方法計(jì)算梢渦流場(chǎng)越來越成為如今的主流方法。但值得注意的是，梢渦流場(chǎng)的速度梯度大，流場(chǎng)分布復(fù)雜。選取合適的湍流模型和對(duì)梢渦流場(chǎng)的網(wǎng)格進(jìn)行細(xì)化處理是十分必要的?，F(xiàn)在廣泛運(yùn)用的RANS方法對(duì)雷諾項(xiàng)進(jìn)行了平均處理，只能得到時(shí)均化后的湍流結(jié)果，忽略了小尺度的湍流結(jié)構(gòu)變化，使得RANS在模擬梢渦中效果很不理想［5］。在計(jì)算網(wǎng)格方面，Turnock等［6］提出了通過尋找壓力最低點(diǎn)定位梢渦軌跡的方法，通過該方法極大地提高了計(jì)算精度。Zhang等［7］則認(rèn)為數(shù)值計(jì)算中渦核徑向網(wǎng)格節(jié)點(diǎn)數(shù)應(yīng)大于15才能滿足計(jì)算精度。

本文的計(jì)算模型是剖面為NACA 16020翼形的半橢圓形水翼，使用不同的湍流模型計(jì)算其流場(chǎng)分布，并與實(shí)驗(yàn)值做對(duì)比，以選取合適的湍流模型。在此基礎(chǔ)上，使用氣泡靜力方程計(jì)算初始梢渦空泡數(shù)。

1 數(shù)值方法

1.1 RANS模型介紹

RANS模型是一種將流體的質(zhì)量、動(dòng)量等時(shí)均化處理的湍流模型［8］。雷諾平均方法不用求解各個(gè)尺度的湍流流動(dòng)，只計(jì)算平均運(yùn)動(dòng)。當(dāng)流體的湍流結(jié)構(gòu)較為復(fù)雜時(shí)，RANS模型只能給出時(shí)均化后的結(jié)果。所以，對(duì)于復(fù)雜的流動(dòng)現(xiàn)象，RANS模型是不適用的。

雷諾平均后的N-S方程如下所示：

式中：xi(i=1,2,3)，xj(j=1,2,3)，xl(l=1,2,3)為三維笛卡爾坐標(biāo)系下的方向坐標(biāo)；t為時(shí)間；ui，uj和ul為速度分量；ρ為流體密度；u為動(dòng)力粘性系數(shù)；為雷諾應(yīng)力；fi為體積力；

1.2 LES模型介紹

湍流流動(dòng)中包含著各種尺度的湍流結(jié)構(gòu)，大尺度渦主要指尺寸大于平均流動(dòng)（注：剪切層厚度）的湍流結(jié)構(gòu)。與大尺度渦相比，小尺度渦主要起耗散湍流能量的作用?；谠摶粳F(xiàn)象，LES使用直接數(shù)值求解的方法計(jì)算大尺度渦，建立模型求解小尺度渦。模型中分離大、小尺度渦的分界尺度稱為過濾尺度，用Δ表示。它相較于普通的RANS模型，要求更細(xì)致的網(wǎng)格分布和更多的計(jì)算資源［9］。

LES的控制方程如下所示：

式中：Lij為L(zhǎng)eonard應(yīng)力，代表著大尺度渦之間的相互作用；Cij為交互應(yīng)力，代表著大尺度渦和小尺度渦之間的相互作用；Rij為亞網(wǎng)格雷諾應(yīng)力，代表著小尺度渦之間的相互作用［10］。Lij，Cij和Rij分別表示如下：

1.3 DES模型介紹

DES模型是一種介于RANS和LES之間的湍流模型，它在邊界層內(nèi)使用RANS模型，其他區(qū)域使用LES模型，所占的計(jì)算資源比RANS模型多，但比LES模型少。本文研究的DES為基于k-ω的湍流模型。

該DES模型與k-ω模型相比，對(duì)耗散項(xiàng)Yk提出了如下修正：

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

本文的計(jì)算模型是剖面為NACA 16020翼形的半橢圓形水翼，實(shí)驗(yàn)［3］已在空泡水洞中完成，來流速度為10 m/s，翼攻角為10°。入口距離水翼3倍弦長(zhǎng)，設(shè)置為速度入口；出口距離水翼10倍弦長(zhǎng)，設(shè)置為壓力出口。其中，原點(diǎn)設(shè)置在水翼根部中心位置。所建計(jì)算域如圖1所示。

圖1 計(jì)算域Fig.1 Computational domain

根據(jù)計(jì)算域的幾何形狀，本文主要用H型網(wǎng)格對(duì)計(jì)算域進(jìn)行整體網(wǎng)格劃分?？紤]到梢渦的形成與翼邊界層內(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ì)量，第1層網(wǎng)格尺度y+在1～30之間。為與實(shí)驗(yàn)結(jié)果相對(duì)比，本文建立的計(jì)算域與實(shí)驗(yàn)一致。本文分別建立了2個(gè)坐標(biāo)系，其中絕對(duì)坐標(biāo)系的原點(diǎn)設(shè)在水翼底部弦長(zhǎng)中心處，x，y，z的方向如圖1所示。為便于對(duì)梢渦進(jìn)行分析，本文還建立了相對(duì)坐標(biāo)系，其原點(diǎn)位于各個(gè)界面處梢渦渦核的中心，x，y，z的方向與絕對(duì)坐標(biāo)系相同。在沒有明確指出時(shí)，默認(rèn)為絕對(duì)坐標(biāo)系。

根據(jù)Turnock等［6］提供的方法：通過初步計(jì)算得到尾流區(qū)域每個(gè)截面壓力最低點(diǎn)的位置，可以認(rèn)為這些點(diǎn)的連線即是梢渦軌跡，并對(duì)該連線進(jìn)行徑向和軸向的網(wǎng)格加密，得到優(yōu)化網(wǎng)格。同時(shí)，Dehghan等［8］也指出，在數(shù)值計(jì)算中，渦核徑向網(wǎng)格節(jié)點(diǎn)數(shù)應(yīng)大于15才能滿足計(jì)算精度。劃分渦核內(nèi)網(wǎng)格節(jié)點(diǎn)為20×20，總網(wǎng)格數(shù)為800萬。

3 梢渦流場(chǎng)數(shù)值研究

三維水翼梢渦處流場(chǎng)結(jié)構(gòu)復(fù)雜，主渦和二次渦的相互干擾影響為數(shù)值模擬帶來了相當(dāng)?shù)碾y度，同時(shí)，梢渦的能量耗散問題也給湍流模型的使用造成了困難。在應(yīng)用不同的湍流模型時(shí)，軸向速度與實(shí)驗(yàn)值相比，存在著較大的不同。Spalart等［11］指出，經(jīng)過旋轉(zhuǎn)和曲率修正的湍流模型能夠有效抑制渦核處的粘性耗散，從而提高計(jì)算精度。但經(jīng)計(jì)算發(fā)現(xiàn)，曲率修正后的湍流模型會(huì)反向減小梢渦處的軸向速度，使其與實(shí)驗(yàn)值相差更大。圖2和圖3給出了x方向3個(gè)位置處梢渦流場(chǎng)無因次軸向速度(U/U∞)和切向速度(Vt/U∞)沿y軸分布的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果。其中，軸向是指與梢渦渦線平行的方向，切向?yàn)榕c梢渦渦線垂直的方向，由于梢渦軌跡近似于x軸方向（來流方向）平行，這里認(rèn)為軸向速度為x方向速度，切向速度為z方向速度。

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

由圖2可知，k-ω模型和基于k-ω的DES模型在預(yù)報(bào)梢渦軸向速度方面存在較大誤差，特別是渦核處的速度分布，與實(shí)驗(yàn)值截然相反。這可能是由于這兩個(gè)模型都過大地預(yù)報(bào)了梢渦處的粘性耗散，使得速度下降得過快。而使用曲率修正后的湍流模型則更加大了粘性耗散，這明顯與文獻(xiàn)［10］的結(jié)論相異，可能的原因是曲率修正的使用建立在一些先決條件上，否則會(huì)得到相反的結(jié)論，對(duì)曲率修正的研究會(huì)在以后的工作中繼續(xù)展開。LES湍流模型則給出了較好的預(yù)報(bào)結(jié)果，特別是x/C=0.1處渦核的速度分布與實(shí)驗(yàn)吻合很好。而在梢渦空泡觀測(cè)中，發(fā)現(xiàn)梢渦空泡的初生位置一般在x/C=0.1附近。所以對(duì)該位置速度流場(chǎng)的準(zhǔn)確模擬為下一步研究梢渦空泡的初生打下了良好的基礎(chǔ)。

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

圖4 梢渦軌跡方向最小壓力系數(shù)分布Fig.4 Pressure coefficients along the tip vortex trajectory

由圖3可知，各湍流模型對(duì)梢渦切向速度分布的預(yù)報(bào)效果均較好，特別是LES模型在x/C=0.3處梢渦粘性出現(xiàn)較大耗散的情況下也能與實(shí)驗(yàn)值符合較好。但在局部坐標(biāo)z＜0的區(qū)域，各湍流模型預(yù)報(bào)的結(jié)果都與實(shí)驗(yàn)值存在較小的誤差，這可能是由于三維水翼壁面對(duì)主渦的干擾所致。圖4給出了沿x方向渦核處的最低壓力系數(shù)分布（Cp），LES，DES和RANS模型的最小壓力都出現(xiàn)在x/C= 0.355 5處，對(duì)應(yīng)的最小壓力系數(shù)分別為-7.38，-7.05和-5.98。很明顯，RANS湍流模型過高地預(yù)報(bào)了粘性耗散，使最小壓力系數(shù)與LES和DES的預(yù)報(bào)結(jié)果差距較大。從圖4可以明顯看出，LES和DES湍流模型給出的最低壓力系數(shù)曲線存在一個(gè)較為明顯的二次降低區(qū)域。該區(qū)域的存在是由于主渦沿梢渦軌跡方向行進(jìn)時(shí)，會(huì)將水翼產(chǎn)生的其他尾渦卷曲進(jìn)來，當(dāng)兩渦螺旋式纏繞并合并時(shí)，會(huì)增大主渦的渦強(qiáng)，從而使渦核內(nèi)流速增大，壓強(qiáng)降低；與此同時(shí)，在主渦的行進(jìn)過程中又會(huì)不斷地耗散能量，所以存在一個(gè)相對(duì)平衡的點(diǎn)，使主渦的耗散能量和卷曲獲得的能量相等。

4 氣泡靜力平衡和初始梢渦空泡數(shù)計(jì)算

描述微氣泡體積變化的經(jīng)典Rayleigh-Plesset方程為

等式左邊為動(dòng)態(tài)項(xiàng)，代表了氣泡體積的變化；右邊為靜態(tài)項(xiàng)，僅與氣泡的受力情況有關(guān)。當(dāng)氣泡體積變化較為緩慢時(shí)，可忽略方程左邊的動(dòng)態(tài)項(xiàng)，方程簡(jiǎn)化為研究空泡初生時(shí)的靜力平衡方程

式（8）和式（9）中：Pe為氣泡周圍的環(huán)境壓力，當(dāng)氣泡的體積很小時(shí)，Pe通常取氣泡中心點(diǎn)的壓力；Pv為飽和蒸汽壓；Pg為氣泡中的氣體壓力；Rb為氣泡半徑；S為氣泡表面的張力系數(shù)；R0為氣泡初始半徑；v為運(yùn)動(dòng)粘性系數(shù)；pv和pencoutner為氣泡周圍壓力；U為流體速度；Ub為氣泡速度；pg0為初始?xì)馀輧?nèi)氣體壓力。

當(dāng)氣泡周圍的環(huán)境壓力變化時(shí)，氣泡的半徑會(huì)根據(jù)其平衡狀態(tài)隨之改變。由于發(fā)生空化時(shí)液體氣化的速度非?？?，設(shè)定氣泡中蒸汽壓力恒為常數(shù)。另一方面，水中氣體擴(kuò)散的速度相較于空化的速度來說十分緩慢，所以可以假設(shè)氣泡中氣體的質(zhì)量恒為常數(shù)。因此Pg的表達(dá)式為

最終的氣泡平衡方程為

其中，

圖5給出了球型氣泡的靜力平衡曲線［12］，其中Critical values為臨界邊界線，該邊界線以下的區(qū)域?qū)儆诜瞧胶鈪^(qū)域，處于該區(qū)域的氣泡會(huì)迅速長(zhǎng)大。因此，由邊界線與靜力平衡曲線的交點(diǎn)即可確定臨界壓力與臨界半徑。

圖5 氣泡靜力平衡曲線［12］Fig.5 Static equilibrium curves of bubbles［12］

求得的臨界壓力Pc與臨界半徑Rc如下所示：

當(dāng)流場(chǎng)壓力小于臨界壓力Pc時(shí)，氣泡會(huì)迅速長(zhǎng)大，梢渦空泡隨之產(chǎn)生。當(dāng)使用氣泡平衡方程時(shí)，不僅要考慮流場(chǎng)環(huán)境，還需要考慮氣核的初始半徑。這里分別設(shè)定氣泡的初始半徑R0=30，50，100 μm，Pv=3 540 Pa，研究不同初始半徑下的梢渦空泡數(shù)。

已知三維水翼的最小壓力系數(shù)為-7.38（使用LES模型計(jì)算的結(jié)果），當(dāng)梢渦空泡初生時(shí)，Pe=Pc。聯(lián)立式（13）和式（14），可計(jì)算出Rc，Pc和初生空泡數(shù)δi，計(jì)算結(jié)果如表1所示。

表1 三維水翼初始梢渦空泡計(jì)算結(jié)果Table 1 The computational results of hydrofoil vortex cavitation inception

從以上計(jì)算結(jié)果可以看出，不同初始半徑下臨界壓力相差較大，導(dǎo)致各個(gè)不同初始半徑的初始梢渦空泡數(shù)相差很大。從以上結(jié)果還可以發(fā)現(xiàn)，臨界壓力總是小于汽化壓力，所以當(dāng)對(duì)控制梢渦空泡初生要求較高時(shí)，可使用偏于安全的汽化壓力來計(jì)算初始梢渦空泡數(shù)。但在需要對(duì)初始梢渦空泡數(shù)進(jìn)行精確計(jì)算時(shí)，則必須使用臨界壓力來計(jì)算梢渦空泡數(shù)。

使用氣泡靜力方程能準(zhǔn)確計(jì)算出初始梢渦空泡數(shù)，該方法為后續(xù)初始梢渦空泡數(shù)的尺度效應(yīng)研究奠定了基礎(chǔ)。

5 結(jié) 論

本文應(yīng)用k-ω，DES和LES模型分別計(jì)算了翼剖面為NACA 16020翼型的橢圓形水翼梢渦流場(chǎng)分布，介紹了氣泡靜力平衡方程，并結(jié)合氣泡靜力平衡方程研究了該水翼的初始梢渦空泡，得到以下結(jié)論：

1）RANS湍流模型過大地預(yù)報(bào)了梢渦的粘性耗散，在應(yīng)用曲率修正以后，結(jié)果并沒有得到改善，反而反向減小了軸向速度，這與文獻(xiàn)［10］的結(jié)論相悖。使用LES湍流模型則能夠得到較好的結(jié)果。

2）應(yīng)用LES湍流模型能準(zhǔn)確預(yù)報(bào)出由于尾渦的卷曲合并對(duì)梢渦處壓力場(chǎng)的影響，為后續(xù)的梢渦空泡研究奠定了基礎(chǔ)。

3）在使用氣泡靜力平衡方程預(yù)報(bào)初始梢渦空泡數(shù)的情況下，當(dāng)對(duì)控制梢渦空泡初生的要求較高時(shí)，可使用偏于安全的汽化壓力來計(jì)算初始梢渦空泡數(shù)。在需要對(duì)初始梢渦空泡數(shù)進(jìn)行精確計(jì)算時(shí)，必須使用臨界壓力計(jì)算梢渦空泡數(shù)。

［1］ BESANT W H.A treatise on hydrostatics and hydrody?namics［M］.London：Deighton，Bell，1859.

［2］ KELLER A P.Cavitation scale effects-empirically found relations and the correlation of cavitation number and hydrodynamic coefficients［C］//CAV 2001：Fourth International Symposium on Cavitation.Pasadena，CA，USA：California Institute of Technology，2001.

［3］ FRUMAN D H，DUGUé C，PAUCHET A，et al.Tip vortex roll-up and cavitation［C］//Proceedings of the 19th Symposium on Naval Hydrodynamics.Washing?ton D C，USA：National Academy，1992：633-654.

［4］ ARNDT R E A，ARAKERI V H，HIGUCHI H.Some observations of tip-vortex cavitation［J］.Journal of Flu?id Mechanics，1991，229：269-289.

［5］ 吳瓊，葉騫，孟國(guó)香.不同湍流模型在旋流非接觸搬運(yùn)器仿真中的應(yīng)用和對(duì)比［J］.上海交通大學(xué)學(xué)報(bào)，2013，47（3）：423-427. WU Q，YE Q，MENG G X.Application and compari?son of the different turbulent models simulation for non-contact vortex gripper［J］.Journal of Shanghai Ji?aotong University，2013，47（3）：423-427（in Chi?nese）.

［6］ TURNOCK S R，PASHIAS C，ROGERS E.Flow fea?ture identification for capture of propeller tip vortex evolution［C］//Proceedings of the 26th Symposium on Naval Hydrodynamics.Rome，Italy：Italian Ship Mod?el Basin，Office of Naval Research，2006.

［7］ ZHANG L X，ZHANG N，PENG X X，et al.A review of studies of mechanism and prediction of tip vortex cavitation inception［J］.Journal of Hydrodynamics：Ser.B，2015，27（4）：488-495.

［8］ DEHGHAN M，TABRIZI H B.Turbulence effects on the granular model of particle motion in a boundary lay?er flow［J］.The Canadian Journal of Chemical Engi?neering，2014，92（1）：189-195.

［9］ 甘文彪，周洲，許曉平，等.基于改進(jìn)SST模型的分離流動(dòng)數(shù)值模擬［J］.推進(jìn)技術(shù)，2013，34（5）：595-602. GAN W B，ZHOU Z，XU X P，et al.Investigation on improving the capability of predicting separation in modified SST turbulence model［J］.Journal of Propul?sion Technology，2013，34（5）：595-602（in Chi?nese）.

［10］ MORGUT M，NOBILE E，BILU? I.Comparison of mass transfer models for the numerical prediction of sheet cavitation around a hydrofoil［J］.International Journal of Multiphase Flow，2011，37（6）：620-626.

［11］ SPALART P R，SHUR M.On the sensitization of tur?bulence models to rotation and curvature［J］.Aero?space Science and Technology，1997，5（1）：297-302.

［12］ CHAHINE G L.Nuclei effects on cavitation inception and noise［C］//Proceedings of the 25th Symposium on Naval Hydrodynamics.St.John's，Newfoundland and Labrador，Canada：［s.n.］，2004.

Numerical study of hydrofoil tip vortex fluid field

PU Jijun，XIONG Ying
Department of Naval Architecture Engineering，Naval University of Engineering，Wuhan 430033，China

Three different models，k-ω，DES and LES，are conducted in the analysis of the tip vortex flow field.In order to reduce the discrete error induced by the grid，mesh refinement is applied to the area of the tip vortex core in numerical simulations.The axis and tangential velocities of the tip vortex flow field with no cavitation are calculated，and the calculated velocities agree well with the experimental results.On the basis of this process，the influence of vortex roll-up on the tip vortex pressure filed is discussed，and bubble static equilibrium is proposed by which the tip vortex cavitation inception number is computed.

tip vortex；bubble static equilibrium；cavitation inception；hydrofoil

U661.313

A

10.3969/j.issn.1673-3185.2017.01.002

2016-06-03

2016-12-28 16:05

國(guó)家自然科學(xué)基金資助項(xiàng)目（51179198）

蒲汲君，男，1991年生，碩士生。研究方向：流體力學(xué)。E-mail：211982361@qq.com熊鷹（通信作者），男，1958年生，博士，教授，博士生導(dǎo)師。研究方向：船舶流體力學(xué)。E-mail：ying_xiong28@126.com

http://www.cnki.net/kcms/detail/42.1755.TJ.20161228.1605.038.html期刊網(wǎng)址：www.ship-research.com

蒲汲君，熊鷹.三維水翼梢渦流場(chǎng)數(shù)值研究［J］.中國(guó)艦船研究，2017，12（1）：8-13，26. PU J J，XIONG Y.Numerical study of hydrofoil tip vortex fluid field［J］.Chinese Journal of Ship Research，2017，12（1）：8-13，26.