螺旋槳梢渦空泡初生及尺度效應研究

熊 鷹，韓寶玉，時立攀

（海軍工程大學 船舶與動力學院，武漢430033）

1 引 言

近年來，氣泡動力學理論結合粘性流數值方法預報梢渦空泡初生問題得到迅速發展。Hsiao等[1-3]對該理論進行了完善，考慮了氣泡與流體之間的滑移速度影響，并應用到螺旋槳梢渦空泡初生預報中。計算過程中，作者用氣泡表面平均壓力代替氣泡中心點處壓力，提出了SAP球形氣泡模型，得到了令人滿意的結果。應用氣泡動力學理論預報螺旋槳梢渦空泡初生的關鍵問題是對螺旋槳梢渦流場的數值模擬。

應用RANS方法計算梢渦流場的難點是數值離散精度及湍流模型。Spall[4]通過對網格的特殊處理，認為在梢渦流場計算中只要渦核處網格尺度合適，二階離散的精度就可以達到工程要求。在湍流模型研究方面，韓寶玉等[5]通過對兩方程和代數雷諾應力模型（EARSM）進行旋轉和曲率修正，研究了湍流模型對機翼梢渦流場數值模擬的影響。發現經過旋轉和曲率修正的EARSM湍流模型計算結果與試驗吻合較好，而未經修正的湍流模型則過高地預報了渦核處的粘性耗散。

應用RANS方法準確計算實尺度螺旋槳梢渦流場時，由于網格數較大，實用性受到限制。而應用模型尺度螺旋槳計算結果換算到實尺度螺旋槳需要考慮尺度效應的影響?？张莩跎c流場中壓力分布密切相關，通常認為當流場中最小壓力達到或低于某臨界值時空泡初生，因此在計算空泡初生空泡數前必須得到流場中最小壓力分布。在工程應用中，通常假定螺旋槳空泡初生空泡數σi與最小壓力系數Cpmin的關系式為σi=-Cpmin，前人對空化初生尺度效應的研究大多數基于上述假定。

根據上述假設及McCormick建立的不同尺度流場最小壓力系數關系式，可得：

Shen等[6]應用線渦理論推導了梢渦中最小壓力系數的表達式，得出n的值為0.4。Bulten[7]以DTRC 4119槳為研究對象，根據最小壓力系數得到梢渦空泡初生空泡數，分析了梢渦空泡初生尺度效應問題，發現當進速系數不同時指數n的值不同，但均在0.2左右，作者建議n的值取為0.35，可以保證一定的安全余量。Amoromin[8]系統分析了渦流中尺度效應問題，給出了n的兩種取值情況，建議在層流中n取為0.4，在湍流中n取為0.24。Hsiao等[3]應用球形氣泡表面平均壓力模型研究了螺旋槳空泡初生的尺度效應問題，近尾流梢渦流場應用直接數值模擬方法求解，除出口邊界條件外，其他邊界條件通過RANS方法獲得。在計算氣泡運動時，合理考慮了流場中初始氣泡的分布，并用聲學準則判斷空泡初生。尺度效應計算結果顯示，不同尺度間指數n的取值有一定差別，從小尺度到中等尺度換算時，n的值為0.22，與Amoromin推薦的湍流中n的取值相近，而中等尺度到大尺度換算時n的值較小，作者認為原因是大尺度模型網格劃分不夠精細導致的計算誤差。Edge[9]基于簡化的Reyleigh-Plesset方程推導了適用于高速射流的空泡初生空泡數尺度換算公式，將計算結果與應用氣泡動力學理論計算出的空泡初生空泡數進行了對比，發現當氣核初始尺寸較大時，兩種方法的計算結果吻合較好。

本文建立了氣泡動力學數值模型并對DTMB 5168槳梢渦空化初生進行了預報。螺旋槳梢渦流場應用RANS方法計算，湍流模型為經過旋轉和曲率修正的代數雷諾應力模型（EARSM-CC）。并從微氣泡動力學角度推導了不同尺度模型空化初生空泡數換算公式，建立了空化初生尺度換算模型，應用不同尺度螺旋槳模型計算結果對換算模型相關參數進行了確定，對螺旋槳空化初生尺度效應進行了研究。

2 計算方法

2.1 數值方法

控制方程為雷諾平均N-S方程，計算中認為流體不可壓縮。對方程的離散采用有限體積法，離散精度為二階，采取的算法為全隱式多網格耦合算法。實驗表明梢渦系統的高旋度及流線的高曲率造成的近似剛體旋轉運動穩定了梢渦內流場分布，使梢渦尾流中的湍動能迅速衰減。為了捕捉渦的這種穩定效應，湍流模型采用代數雷諾應力模型[10]（EARSM），計算中考慮了Wallin等[11]提出的旋轉和曲率修正方法。

2.2 氣泡運動模型

球形氣泡在流場中的運動包括體積脈動和位移變化兩種形式。為了準確模擬球形氣泡在壓力場的體積脈動，Hsiao等[1]對經典的Reyleigh-Plesset方程進行了修正，考慮了氣泡與流體之間的滑移速度影響，得到如下控制方程：

式中：rb為氣泡半徑，pv為流體飽和蒸汽壓，k為空氣多變指數，r0為氣泡的初始半徑，pe為氣泡周圍的環境壓力，S為氣泡表面張力系數，μ為流體動力粘性系數，和分別為流場和氣泡的速度。氣泡的初始壓力pg0可由下式得出：

在求解方程（2）時，經典的氣泡模型認為pe為氣泡中心點處的壓力。該簡化對于渦核處流場壓力小于飽和蒸汽壓時，所求得的氣泡半徑會無限變大，這與實際情況是不相符的。為此，本文應用氣泡表面受到的平均壓力代替氣泡中心點處的壓力。

為了計算氣泡運動過程中的位移變化，需要引入氣泡運動微分方程式：

式中：阻力系數CD由下式確定：

氣泡脈動在流場距離氣泡中心點l（l?rb）處輻射聲壓可由Fitzpatrick-Strasberg方程確定：

噪聲級SPL可由下式求得：

式中：pref=10-6N/m2

球形氣泡在水翼流場中的運動軌跡和體積脈動通過求解方程（2）和方程（4），求解該方程組的數值方法為四階Runge-Kutta法，方程組中流場的速度和壓力變量應用（2.1）節建立的數值方法求解。應用RANS方法求解流場分布時，只能得到各變量在網格節點處的值，而氣泡的運動位置是不定的，氣泡運動位置處的變量數值可根據網格節點進行插值得到。本文對計算域采用全六面體網格劃分，流場中任意點（x，y，z）處變量值用該點所在的網格八個節點處變量插值得到，其表達式為：

其中

（xi，yi，zi）為點（x，y，z）所在網格的節點坐標，φ，ψ，φ 為插值系數，如果氣泡位置處于某網格中，則求得的 φ，ψ，φ 滿足0≤φ≤1，0≤ψ≤1，0≤φ≤1。求解方程（10）的數值方法為Newton-Raphson法。

數值計算中，方程（10）中的pe為運動氣泡壁面所受的平均壓力，該壓力值可通過高斯積分求得。當氣泡半徑達到一定數值時，氣泡表面點可能超出了壁面，計算過程中，超出點處的環境壓力用距離該點最近的壁面處網格中心點處環境壓力值表示。

2.3 空泡初生尺度換算模型

Reyleigh-Plesset控制方程為非線性二階常微分方程，完全考慮方程各項的影響是非常復雜的。文獻[9]從數值計算的角度探討了經過一定簡化的R-P方程與完整R-P方程計算結果的差異，發現在低頻壓力場中，對R-P方程左側動態項作近似零處理對計算結果的影響很小。螺旋槳梢渦內的壓力場屬低頻壓力場，因此在研究中可以對R-P方程進行簡化。其簡化形式為：

假定當氣泡周圍平均環境壓力pe達到臨界值pc時，空泡半徑rb達到可見尺寸rv。從（12）式可以看出，在同一流體中，當空泡達到可見尺寸時，對不同尺度模型而言方程右側保持不變，即pc-pv保持不變。其表達式為：

空泡數σi可表示為：

空化初生時，臨界壓力pc等于流場最小壓力pmin，可得：

選定兩個不同尺度模型梢渦流場進行分析，定義兩個尺度為模型尺度和實尺度，梢渦空化初生空泡數σim和σis可表示為：

式中：下標m和s分別表示模型尺度和實尺度，由于pc-pv為常數，將（16）式代入（17）式可得：

模型尺度和實尺度梢渦流場最小壓力系數的關系式為：

將（19）式代入（18）式可得梢渦空泡初生尺度換算公式：

2.4 計算域及網格劃分

計算對象為均勻流中的螺旋槳模型，槳模為泰勒水池DTMB 5168五葉模型槳。為與實驗[12]一致，設定槳模直徑D=0.402 7m。由于槳轂形狀對梢渦流場的影響較小，本文用貫穿計算域的圓柱面對槳轂進行簡化，有助于形成較高質量的網格。選取單個槳葉所在的圓柱單通道作為計算域，計算域的周期性對稱面位于兩槳葉中間，進口與槳盤面距離為2D，出口與槳盤面距離為5D，外圓柱面直徑為2.5D。如圖1所示。

圖1 計算域Fig.1 The calculation domain

采用全六面體網格形式對螺旋槳計算域進行網格劃分，整個計算域可看做一個大塊，整體的網格為H-H型。為了準確模擬螺旋槳葉片邊界層及其附近流場內流動情況，對槳葉附近切割出的小塊用O型網格進行處理，網格總數為443萬。螺旋槳葉片處第一層網格尺度設定為2.5×10-5D，保證葉片處y+在1-10之間。由于梢渦渦核內徑向速度梯度較大，為了減小該區域的數值離散誤差，網格生成時特意對梢渦渦核區域的網格進行加密，保證渦核內每個方向上網格節點數達到20以上。計算域及部分區域網格見圖2。

在邊界條件設置上，由于計算涉及到的進速系數J=U0/nD=1.1（U0為軸向來流速度），對應的螺旋槳雷諾數Rn=nD2/ν=3.92×106。設定進口和外圓柱面為速度進口，方向為x軸向，其大小與旋轉坐標系的轉速設定值（旋轉軸為x軸，旋轉方向滿足右手法則）滿足上述雷諾數，入口處的湍流強度為5%，渦粘比為10；出口設為壓力邊界條件，平均靜壓為1atm。螺旋槳葉表面為旋轉、不可滑移物面邊界條件，旋轉速度與旋轉坐標系相同。計算域兩側面設定為旋轉周期性邊界條件。

圖2 計算網格Fig.2 The calculation mesh

3 計算結果及其分析

3.1 尾流場計算結果

對螺旋槳尾流場的數值模擬不僅有助于理解梢渦的發生和發展情況，更由于尾流對槳后舵及其他附體的干擾，使其成為設計這些裝置的關鍵因素。同時螺旋槳尾流場計算對于分析預報螺旋槳激振力至關重要。圖3為進速系數J=1.1時，螺旋槳下游x/R=0.238 6處部分區域三個速度分量云圖（軸向Vx，切向Vt，徑向Vr）的計算結果與實驗結果，其中實驗結果不包含槳轂附近的流場分布云圖。比較發現，在梢渦區域，槳葉后尾流場及槳葉間流場上的三個速度分量計算結果均與實驗結果吻合，為后續梢渦流場相關研究提供了理論基礎和支撐。

在尾流場中，梢渦流場是本文研究的重點。為了驗證數值方法對螺旋槳梢渦預報的準確度，本文對螺旋槳下游x/R=0.238 6，渦核軸線位置所處半徑r處一定角度θ范圍內的速度Vx、Vt、Vr進行了計算。根據實驗結果渦核軸線位置半徑r/R=0.92，計算結果與實驗值的比較見圖4。其中渦核軸線處周向角θ定為0。

圖3 x/R=0.238 6處剖面速度云圖Fig.3 Velocity contour on cross plane at x/R=0.238 6

圖4中，切向速度分布曲線從左至右第一個波谷是葉片后尾流位置，第二個波谷是梢渦渦核位置?？梢钥闯霰疚牡挠嬎阒蹬c實驗結果是吻合的，這說明計算中網格在渦核位置y方向×z方向節點數達到20×20是足夠達到精度的，同時經過旋轉和曲率修正的EARSM湍流模型有效解決了常用數值方法中渦核內粘性耗散計算值過大的問題。

3.2 梢渦空泡初生計算

根據進速系數J=1.1時的螺旋槳梢渦流場分布計算結果，通過在適當位置釋放初始半徑r0=10 μm，50 μm，100 μm的氣泡，計算了在不同空泡數流場中不同初始尺寸氣泡的運動形態。計算中，氣核釋放位置x方向在槳盤面前方距離△x=0.1R，而y方向和z方向釋放位置需要反復試驗，以保證氣核在最小壓力點前進入渦核。為了保證氣核運動軌跡的準確性及計算的穩定性，應選擇足夠小的時間步，計算中時間步選為1×10-8s。

在計算氣泡的運動形態時，方程（2）中pe為運動氣泡壁面所受的平均壓力，由于計算選取的螺旋槳經過葉梢卸載處理，梢渦渦核尺寸較小，當氣泡直徑達到初生判定標準時，平均壓力值與梢渦的最小壓力差別較大，導致螺旋槳空化初生空泡數計算值可能遠小于-Cpmin。根據水翼的研究結果，在最小壓力系數Cpmin=-3.47，空化初生空泡數σi=3.4時，梢渦流場中壓力p=233 7 Pa的等值面橫截面最大尺寸在0.75 mm左右。以此為依據，選取σ=3.1，3.0和2.9三個空泡數時的流場，分別對梢渦空泡初生過程中不同尺度微氣泡的運動形態進行了計算分析，計算中取距離氣泡中心點1 m處輻射聲壓值為研究對象。圖5為σ=2.9，r0=100 μm時氣泡運動過程中氣泡半徑，氣泡中心點處流場壓力及輻射聲壓隨時間變化圖。根據計算結果可以從光學和聲學角度判斷空化是否初生，應用光學準則時，根據肉眼的分辨率，認為氣泡直徑達到1 mm時空泡初生，應用聲學準則時，本文假定當距離氣泡中心點1 m處輻射聲壓達到10 Pa時能檢測到聲壓。

圖4 x/R=0.238 6處通過渦核軸線切向速度分布Fig.4 Velocity distribution across the vortex core axis in the tangential direction at x/R=0.238 6

圖5 σ=2.9，r0=100 μm時氣泡的運動Fig.5 Sub-visual cavitation at σ=2.9 and r0=100 μm

計算結果見表1和表2，從計算結果可以看出，根據光學準則，在空泡初生前，r0越大的氣泡在運動至渦核低壓區時氣泡半徑越大，當空泡初生時，r0越小的氣泡運動至渦核低壓區時氣泡半徑越大，其原因是由于渦核尺寸較小，計算中r0較小的氣泡在渦核內所受的平均壓力大于r0大的氣泡，導致r0較小的氣泡壁面加速度較大。

根據光學準則和聲學準則，螺旋槳梢渦空泡初生空泡數σi計算結果分別為3.0和2.9，高于實驗觀測結果2.0。原因是本文計算中采用的球形氣泡模型只適用于空泡初生過程，且能夠觀測到空泡初生的時間很短，而實驗中能夠看到相對而言形態比較穩定的空泡流時才認為梢渦空泡初生，此時空泡數已經低于空泡初生理論空泡數?？张莩跎鷷r，實驗中看到的空泡流一般在葉梢后一定距離處，這可能是大量的氣泡經過低壓區潰滅后在梢渦尾流中匯集而成。

表1 氣泡最大半徑Tab.1 Maximum bubble size

表2 輻射噪聲水平Tab.2 Acoustic noise level

3.3 梢渦空泡初生換算

應用梢渦空泡初生尺度換算公式從模型尺度到實尺度空泡初生空泡數換算時，指數n對計算結果的影響是不可忽略的。為了確定指數n的取值范圍，本文分別選取了三個不同尺度的螺旋槳。計算模型的直徑D，來流速度V0，螺旋槳轉速n，旋轉雷諾數Rn及最小壓力系數的計算結果見表3。

根據水翼和螺旋槳梢渦最小壓力系數計算結果，根據（19）式可得指數n的值為0.12。由于n值越大，計算結果越安全和保守，為了通過尺度換算得到大尺度模型安全保守的空泡初生空泡數，本文建議n的值取為0.2。

首先應用本文介紹的梢渦空泡初生預報方法預報螺旋槳尺度1梢渦空泡初生空泡數，空泡初生判斷準則為光學準則，空泡初生時氣泡直徑為1 mm，然后應用方程（20）預報螺旋槳尺度2和尺度3的空泡初生空泡數，并與應用微氣泡理論計算出的結果進行對比。計算結果見表4，表中應用方程（20）式計算結果和應用微氣泡理論計算結果簡寫為“公式”和“數值方法”。表中應用公式求解時，n的值取為0.2。

從計算結果看，應用本文建立的空化初生尺度換算模型得到的尺度2和尺度3模型梢渦空泡初生空泡數稍高于應用數值方法預報結果，但相差不大。這說明換算結果與數值模擬結果相比是更安全的。

表3 三個尺度螺旋槳相關特征Tab.3 Characteristics of the three propeller scales

表4 螺旋槳空泡初生空泡數預報結果Tab.4 Calculation results of incipient cavitation number of propeller

4 結 論

本文通過求解氣泡動力學方程研究了氣泡在螺旋槳梢渦流場中的運動形態，建立了螺旋槳梢渦空泡初生預報方法，并根據簡化Reyleigh-Plesset方程推導了梢渦空泡初生尺度換算方程，得出如下結論。

（1）初生空泡數數值計算結果高于試驗觀測結果，其原因可能是應用球形氣泡模型判斷空泡初生時，能夠觀測到空泡的時間很短，而實驗中能夠看到相對而言形態比較穩定的空泡流時才認為梢渦空泡初生，此時空泡數已經明顯低于空泡初生理論空泡數。

（2）在空泡初生前，r0越大的氣泡在運動至渦核低壓區時氣泡半徑越大，當空泡初生時，r0越小的氣泡運動至渦核低壓區時氣泡半徑越大，其原因是由于渦核尺寸較小，計算中r0較小的氣泡在渦核內所受的平均壓力大于r0大的氣泡，導致r0較小的氣泡壁面加速度較大。

（3）應用本文建立的空化初生尺度換算模型得到的結果與數值模擬結果基本吻合，而且其結果偏于安全。

[1]Hsiao C T,Pauley L L.Study of tip vortex cavitation inception using Navier-Stokes computation and bubble dynamics model[J].Journal of Fluids Engineering,1999,121:198-204.

[2]Hsiao C T,Jain A,Chahine G L.Effect of gas diffusion on bubble entrainment and dynamics around a propeller[C]//26th Symposium on Naval Hydrodynamic.Washington D C,National Academy Press,2006.

[3]Hsiao C T,Chahine G L.Scaling of tip vortex cavitation inception for a marine open propeller[C].27th Symposium on Naval Hydrodynamic.Washington D C,National Academy Press,2008.

[4]Spall R E.Numerical study of a wing-tip vortex using the Euler equations[J].Journal of Aircraft,2001,38(1):22-27.

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

[6]Shen Y,Chahine G,Hsiao C T．Effects of model size and free stream nuclei on tip vortex cavitation inception[C]//Fourth International Symposium on Cavitation.Wageningen,2001.

[7]Bulten N,Oprea A I.Evaluation of McCormick’s rule for propeller tip cavitation inception based on CFD results[C]//Sixth International Symposium on Cavitation.Wageningen,The Netherlands,2006.

[8]Amromin E.Two-range scaling for vortex cavitation inception[J].Ocean Engineering,2006,33:530-534.

[9]Edge B A.On analytical and numerical study of cavitation scale effects in high-Reynolds number circular jet flows[D].Ph.D.Thesis,The Pennsylvania State University,2007.

[10]Hellsten A.New advanced k-ω turbulence model for high-lift aerodynamics[J].AIAA Journal,2005,43(9):1857-1869.

[11]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.

[12]Chesnakas C,Jessup S.Experimental characterization of propeller tip flow[C].Twenty-Second Symposium on Naval Hydrodynamic.Washington D C:National Academy Press,1998:156-170.