大側斜螺旋槳敞水與空化性能數值模擬

張光明，章 杰

（中國船舶重工集團有限公司第七一〇研究所，湖北 宜昌 443003）

0 引言

近十幾年，隨著計算機技術的快速發展，計算效率得到大幅度提升，加之國內外大量學者的研究工作，基于數值計算的螺旋槳敞水性能預報已經完全達到工程應用精度要求，基本可以替代螺旋槳模型試驗。螺旋槳的空化研究在最近幾年也得到了廣泛的關注，綜合國內外大量文獻及相關的研究可以發現，目前對螺旋槳的空化數值模擬大多基于RANS方法，采用兩方程湍流模型以及基于氣泡動力學方程的兩相流空化模型來完成。Kawamura 等（2005）[1]運用雷諾應力湍流模型和全空泡模型對東京大學開發的某一螺旋槳進行了數值模擬，空泡的位置與形態與公布的試驗結果相吻合；趙鵬飛（2011）[2]采用3種湍流模型即kε-、SST和雷諾應力模型結合 Z-G-B空泡模型對某一大側斜螺旋槳進行計算，得到的推力和扭矩系數在主要工作工況下的相對誤差均在10%以內，其中雷諾應力模型的誤差最小；Ji B等（2011）[3]用CFX軟件采用雷諾應力和SST湍流模型結合Zwart空泡模型，對標SEIUN-MARU大側斜螺旋槳進行了空泡模擬，計算出的由空泡引起的壓力脈動與試驗結果的誤差在20%以內，可以滿足目前工程應用的要求。

在2011年召開的SMP2011會議上，專門安排了螺旋槳和水翼空泡的數值模擬的研討會（Workshop），相關的實驗與CFD計算均選擇PPTC（Postdam Propeller Test Case）槳作為算例槳，該槳屬于大側斜螺旋槳，既有縱斜，又有側斜。相關的實驗由 SVA完成，并在會后公布了相關的測試結果。本文選取PPTC大側斜螺旋槳作為研究對象，采用基于RANS方程的求解方法來計算PPTC螺旋槳的敞水性能和典型狀態下的空泡性能。

1 理論基礎

本文首先采用 RANS方法來計算螺旋槳的敞水性能，應用基于兩相流的RANS方法來計算螺旋槳的空化性能。在計算域內，不可壓縮的牛頓流體都滿足連續性方程和動量方程[4]。用RANS方程描述可以表示為如下。

兩相流下的RANS方程：

式（1）和（2）就是常說的雷諾平均N-S方程，其中，τij為平均粘性剪切力，方程（2）中出現了雷諾應力項，因此需要引入湍流模型來使方程封閉。通常的做法就是引入湍流模型來解決基于RANS方程的湍流流動問題[5]。

在基于RANS方程的方法中，最常用的湍流模型主要有 2種：渦粘模型和雷諾應力模型，而在FLUENT中，Spalart-Allmaras模型、kε-模型和kω-模型都是基于渦粘模型，根據Boussinesq提出的假設，認為雷諾應力與平均速度梯度成正比。

盡管由于各向同性的假設，渦粘模型在很多情況下并非嚴格有效，但是卻已經為很多實際流動問題的計算提供了十分不錯的結果。因此本文的計算也將采用基于渦粘模型的兩方程湍流模型。

在計算螺旋槳的空化性能時，需要應用FLUENT中嵌入的空泡模型。Singhal完整空泡模型、Schneer&Sauer空泡模型和Z-G-B空泡模型都是基于質量運輸方程的空泡模型，其實質就是在質量運輸方程中以水蒸氣的體積分數的變化來表征水和水蒸氣之間的相變[6]。

基本的氣相運輸方程可以表示為

式中：Vρ為氣相的密度；α代表氣相的體積分數；為氣相的速度矢量；Re和RC表示由于空泡生長和潰滅引起的質量變化值。

在空泡動力學里，氣泡的生長和潰滅取決于Rayleigh-Plesset方程，它是一個描述流體中空泡運動的常微分方程，可表示為如下：

式中：RB表示空泡的直徑大小；ρl和νl分別表示液相的密度和運動粘度；pB和p分別表示空泡的表面壓力和遠場壓力值；γ為液相的表面張力。

假設表面張力、粘度以及慣性效應在表面壓力和遠場壓力差值較大時可以忽略，則 Rayleigh-Plesset方程可以簡化為

2 計算模型及計算域

2.1 計算模型

本文的計算模型為PPTC大側斜螺旋槳。具體的幾何參數如表1。

根據PPTC的幾何特征參數以及各葉切面處的型值，通過三維坐標變換公式，得到各個葉切面在三維坐標中坐標值，然后在三維建模軟件中得到其三維圖形。本文采用的是將三維坐標按一定的編排順序，將坐標值導入ICEM軟件，最終得到的PPTC螺旋槳模型如圖1。

表1 PPTC螺旋槳的幾何特征參數Table 1 Dimensions of the PPTC model propeller

圖1 PPTC螺旋槳三維模型Fig.1 Three-dimensional model of PPTC propeller

2.2 計算域及網格劃分

整個計算域為與槳同軸的圓柱形區域，取D為螺旋槳直徑。大內圓柱用于設置局部加密和旋轉坐標系。如圖2所示。

圖2 計算域示意圖Fig.2 Schematic diagram of computational domain

邊界條件的設置如下：

1）入口條件。入口為速度入口，通過改變來流速度大小來得到不同的進速系數。

2）出口條件。出口為自由流出。

3）壁面條件。除螺旋槳的葉面、葉背和葉稍為無滑移壁面，其他均為滑移壁面。

4）Interface。包裹螺旋槳的小圓柱旋轉來模擬螺旋槳的旋轉運動，大圓柱固定不動，交界面處均設為interface。

在數值計算之前要先進行計算域的離散化處理，也就是常說的網格劃分。目前比較常用的是結構網格與非結構網格。

結構網格具有貼體的特性，能夠生成高質量的邊界層，網格區域光滑，計算效率高，對于流場的細節捕捉到位。但是對于幾何外形特別復雜的模型（如螺旋槳）若要劃分結構網格，則要花費大量時間與腦力進行分塊拓撲，同時生成的網格質量也不高，還容易產生負體積，因此結構網格在處理復雜幾何模型的流動問題上并不適用。

非結構網格能夠很好地填充計算域，網格能夠自適應極其復雜的幾何外形，但網格數量一般較大。本文的研究將采用混合網格的劃分模式，結合兩者的優點：包裹著螺旋槳的小圓柱區域內采用非結構網格劃分，以便自適應螺旋槳復雜的外形；小圓柱與大圓柱之間的區域比較規則，則采用結構網格劃分，以提高計算效率；兩者的交界面處采用interface傳遞數據。

圖3 螺旋槳槳葉網格與槳盤面網格分布（B）Fig.3 Grid distribution of blade and disk of propeller（B）

最終，通過上述網格劃分方法，并設置槳葉附近的邊界層，得到3套不同密度的網格：A套網格總數約為 230萬，其中旋轉區域非結構網格約為130萬；B套網格總數約為330萬，其中旋轉區域非結構網格約為200萬；C套網格總數約為410萬，其中旋轉區域非結構網格約為280萬。

3 螺旋槳敞水數值模擬

3.1 不同湍流模型下敞水計算比較

本節根據之前的網格劃分，選取 A套網格進行敞水性能計算。螺旋槳轉速設為 600 r/min，通過改變來流速度得到不同的進速系數。本節設置PPTC的進速系數在0.6～1.4之間變化。湍流模型為常用的5個兩方程湍流模型：標準k-ε模型、RNG k-ε模型、Realizable k-ε模型、標準k-ω模型和SST k-ω模型。

通過相關螺旋槳性能計算理論[7]，對計算結果的分析與處理，我們得到了不同湍流模型下的推力系數KT、扭矩系數10KQ以及敞水效率η，再與公布的試驗數據對比，結果如表2。

表2 不同湍流模型的敞水計算結果Table 2 Open-water calculation results of different turbulence models

將計算結果繪制成曲線圖，將更直觀地與試驗結果進行對比。通過表2、圖4-6可以很明顯地看出，5種湍流模型都能比較準確地計算出PPTC螺旋槳的敞水性能，計算結果已經完全能滿足工程應用。

圖4 推力系數曲線Fig.4 Trust coefficient curve

圖5 扭矩系數曲線（10KQ）Fig.5 Torque coefficient curve（10KQ）

綜合比較后可以看出采用 RNG k-ε湍流的模型計算的結果誤差是最小的。與試驗結果對比，采用 RNG k-ε計算得到的 KT的誤差絕對值小于3.44%，KQ的誤差絕對值小于6.24%，敞水效率的誤差絕對值小于5.57%。其中，在進速系數較大時誤差通常表現的比較大。這是因為在高進速系數時，推力與轉矩的絕對值較小，因而誤差相對低進速比較大也是合理的。因此，本文后面的計算將采用RNG k-ε湍流模型。

圖6 敞水效率曲線Fig.6 Open-water efficiency curve

3.2 網格無關性驗證

依據上一節的螺旋槳在不同湍流模型下的敞水計算結果，可知應用RNG k-ε湍流模型能相對比較準確地計算出螺旋槳的敞水性能參數。本節將采用RNG k-ε湍流模型和A、B、C共3套不同網格數量的計算域來計算螺旋槳敞水性能，通過對比結果來驗證網格的無關性。

計算后結果分析如表3。

將計算敞水結果繪制成曲線圖，將更直觀地驗證網格的無關性：通過對表中數據進行仔細分析比較，可以看出 3套網格對 PPTC螺旋槳的敞水預報結果與試驗結果十分的吻合。只有在進速系數較大時，誤差略微有所變大，但仍然是一個很小的值；同時可以發現計算結果有隨網格數量增大和變的更為精確的趨勢，但不是很明顯。所以可以得出計算結果在網格數量達到A之后，對計算結果的影響就不再明顯。如圖 7。

表3 螺旋槳敞水性能網格因素比較Table 3 Comparison of grid factors of open-water performance

圖7 螺旋槳敞水計算網格無關性驗證Fig.7 Grid independence verification of open-water calculation

4 螺旋槳空化性能計算

根據前節螺旋槳敞水計算的結果，為了在計算空化流動時達到一定的精度，同時又能合理地利用計算資源，此處選擇B網格。選用RNG k- ε湍流模型，結合使用最多Schnerr and Sauer空泡模型來對PPTC大側斜螺旋槳進行空化計算。根據公布的試驗結果，選定的計算狀態為：進速系數J=1.019，空化數nσ=2.024，轉速n=25 r/s。

根據空化數的定義公式：

最終計算得到的邊界條件為：速度入口VA=6.37 m/s，壓力出口 P0=42 381.168 Pa，其它邊界條件同敞水計算。

在上述邊界條件下，分別計算了有空泡和無空泡的結果，結果匯總如下：

在fluent后處理時，建立氣相體積iso等勢面，取氣體體積分數為 0.2，并與其它文獻的計算結果進行對比[8]，得到如圖8所示。

比對計算結果與試驗結果，可以發現片空泡發生的位置與試驗基本相同，尤其是葉根和葉稍處都出現了空化現象；而對于計算的結果的導邊處的空泡，查閱大量文獻發現，對PPTC螺旋槳的計算結果都不同程度的出現了該現象。

將計算的空化狀態下與無空化狀態下的葉面葉背壓力云圖匯總如圖9。

圖8 空泡體積等勢面圖（體積分數0.2）Fig.8 Volume equipotential surface diagrams of cavitation（volume fraction is 0.2）

圖9 葉面葉背壓力分布云圖Fig.9 Pressure distribution nephogram on surface and back of blade

將葉片各個典型半徑截面處的壓力系數繪制成xy-plot圖。

通過對比螺旋槳葉片典型半徑在有空化和無空化工況下的壓力系數曲線圖（圖 10），可以看出壓力系數曲線在葉片上的變化始終時連續的。在0.7r和0.9r處大約有10%的弦長出現了片空化，從0.95r至葉尖則出現了穩定的片空化，這與試驗結果也是相吻合的。

圖10 葉片典型半徑處壓力系數圖Fig.10 External pressure coefficient diagrams of typical radius of blade

5 結束語

本文基于RANS方程，使用FLUENT對PPTC大側斜螺旋槳在均勻來流中的敞水和空化性能進行數值模擬，計算的結果與試驗結果基本吻合，敞水計算結果與試驗結果的誤差總體小于4%，完全滿足工程應用的要求。在空化計算時，本文選取了SVA公布的第一個空泡試驗結果的典型狀態case2.3.1（J=1.019，nσ=2.024），計算結果出現了穩定的片狀空泡，出現的位置與試驗的結果基本吻合，尤其是在葉根和葉稍處。而對于試驗結果觀察到的稍渦空泡，計算結果始終沒有出現，主要原因可能是在稍渦出現的位置的網格不夠密造成的。