三維水翼梢渦流場數值研究

蒲汲君，熊鷹

海軍工程大學艦船工程系，湖北武漢430033

三維水翼梢渦流場數值研究

蒲汲君，熊鷹

海軍工程大學艦船工程系，湖北武漢430033

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

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

0 引 言

自1859年Besant［1］的早期研究以來，梢渦空泡一直是學者們研究的焦點。當一個氣核進入到渦核當中，如果這里的壓力比臨界壓力更小，那么氣核將會迅速長大并形成一個肉眼可見的空泡。空泡在潰滅過程中會釋放強烈的輻射噪聲并造成螺旋槳和船舵的剝蝕。空泡的類型包括梢渦空泡、片空泡、云空泡和游移型空化。大量實驗發現，梢渦渦核中心處的局部壓力最低，除了梢部有卸載處理的三維水翼和螺旋槳外，梢渦空泡在所有類型的空泡中最早出現（Keller［2］在其相關實驗中已有發現），它產生的噪聲會對艦船隱身性產生巨大影響。為抑制梢渦空泡的產生，很有必要準確計算出梢渦流域的流場特性。但是，由于梢渦本身的不穩定性，準確把握梢渦流場的特性還是一個較為困難的問題。解決該問題的有效途徑是進行高精度的數值計算，它能提供梢渦流域流場的具體細節。

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

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

1 數值方法

1.1 RANS模型介紹

RANS模型是一種將流體的質量、動量等時均化處理的湍流模型［8］。雷諾平均方法不用求解各個尺度的湍流流動，只計算平均運動。當流體的湍流結構較為復雜時，RANS模型只能給出時均化后的結果。所以，對于復雜的流動現象，RANS模型是不適用的。

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

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

1.2 LES模型介紹

湍流流動中包含著各種尺度的湍流結構，大尺度渦主要指尺寸大于平均流動（注：剪切層厚度）的湍流結構。與大尺度渦相比，小尺度渦主要起耗散湍流能量的作用。基于該基本現象，LES使用直接數值求解的方法計算大尺度渦，建立模型求解小尺度渦。模型中分離大、小尺度渦的分界尺度稱為過濾尺度，用Δ表示。它相較于普通的RANS模型，要求更細致的網格分布和更多的計算資源［9］。

LES的控制方程如下所示：

式中：Lij為Leonard應力，代表著大尺度渦之間的相互作用；Cij為交互應力，代表著大尺度渦和小尺度渦之間的相互作用；Rij為亞網格雷諾應力，代表著小尺度渦之間的相互作用［10］。Lij，Cij和Rij分別表示如下：

1.3 DES模型介紹

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

該DES模型與k-ω模型相比，對耗散項Yk提出了如下修正：

2 計算模型和網格劃分

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

圖1 計算域Fig.1 Computational domain

根據計算域的幾何形狀，本文主要用H型網格對計算域進行整體網格劃分。考慮到梢渦的形成與翼邊界層內流動有關，為較準確地模擬邊界層內流場，本文應用O型網格對翼壁面附近的網格進行了處理，以提高網格質量，第1層網格尺度y+在1～30之間。為與實驗結果相對比，本文建立的計算域與實驗一致。本文分別建立了2個坐標系，其中絕對坐標系的原點設在水翼底部弦長中心處，x，y，z的方向如圖1所示。為便于對梢渦進行分析，本文還建立了相對坐標系，其原點位于各個界面處梢渦渦核的中心，x，y，z的方向與絕對坐標系相同。在沒有明確指出時，默認為絕對坐標系。

根據Turnock等［6］提供的方法：通過初步計算得到尾流區域每個截面壓力最低點的位置，可以認為這些點的連線即是梢渦軌跡，并對該連線進行徑向和軸向的網格加密，得到優化網格。同時，Dehghan等［8］也指出，在數值計算中，渦核徑向網格節點數應大于15才能滿足計算精度。劃分渦核內網格節點為20×20，總網格數為800萬。

3 梢渦流場數值研究

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

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

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

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

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

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

4 氣泡靜力平衡和初始梢渦空泡數計算

描述微氣泡體積變化的經典Rayleigh-Plesset方程為

等式左邊為動態項，代表了氣泡體積的變化；右邊為靜態項，僅與氣泡的受力情況有關。當氣泡體積變化較為緩慢時，可忽略方程左邊的動態項，方程簡化為研究空泡初生時的靜力平衡方程

式（8）和式（9）中：Pe為氣泡周圍的環境壓力，當氣泡的體積很小時，Pe通常取氣泡中心點的壓力；Pv為飽和蒸汽壓；Pg為氣泡中的氣體壓力；Rb為氣泡半徑；S為氣泡表面的張力系數；R0為氣泡初始半徑；v為運動粘性系數；pv和pencoutner為氣泡周圍壓力；U為流體速度；Ub為氣泡速度；pg0為初始氣泡內氣體壓力。

當氣泡周圍的環境壓力變化時，氣泡的半徑會根據其平衡狀態隨之改變。由于發生空化時液體氣化的速度非常快，設定氣泡中蒸汽壓力恒為常數。另一方面，水中氣體擴散的速度相較于空化的速度來說十分緩慢，所以可以假設氣泡中氣體的質量恒為常數。因此Pg的表達式為

最終的氣泡平衡方程為

其中，

圖5給出了球型氣泡的靜力平衡曲線［12］，其中Critical values為臨界邊界線，該邊界線以下的區域屬于非平衡區域，處于該區域的氣泡會迅速長大。因此，由邊界線與靜力平衡曲線的交點即可確定臨界壓力與臨界半徑。

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

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

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

已知三維水翼的最小壓力系數為-7.38（使用LES模型計算的結果），當梢渦空泡初生時，Pe=Pc。聯立式（13）和式（14），可計算出Rc，Pc和初生空泡數δi，計算結果如表1所示。

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

從以上計算結果可以看出，不同初始半徑下臨界壓力相差較大，導致各個不同初始半徑的初始梢渦空泡數相差很大。從以上結果還可以發現，臨界壓力總是小于汽化壓力，所以當對控制梢渦空泡初生要求較高時，可使用偏于安全的汽化壓力來計算初始梢渦空泡數。但在需要對初始梢渦空泡數進行精確計算時，則必須使用臨界壓力來計算梢渦空泡數。

使用氣泡靜力方程能準確計算出初始梢渦空泡數，該方法為后續初始梢渦空泡數的尺度效應研究奠定了基礎。

5 結 論

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

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

2）應用LES湍流模型能準確預報出由于尾渦的卷曲合并對梢渦處壓力場的影響，為后續的梢渦空泡研究奠定了基礎。

3）在使用氣泡靜力平衡方程預報初始梢渦空泡數的情況下，當對控制梢渦空泡初生的要求較高時，可使用偏于安全的汽化壓力來計算初始梢渦空泡數。在需要對初始梢渦空泡數進行精確計算時，必須使用臨界壓力計算梢渦空泡數。

［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］ 吳瓊，葉騫，孟國香.不同湍流模型在旋流非接觸搬運器仿真中的應用和對比［J］.上海交通大學學報，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］ 甘文彪，周洲，許曉平，等.基于改進SST模型的分離流動數值模擬［J］.推進技術，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

國家自然科學基金資助項目（51179198）

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

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

蒲汲君，熊鷹.三維水翼梢渦流場數值研究［J］.中國艦船研究，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.