螺旋槳梢渦及梢渦空泡數(shù)值模擬

劉芳遠，傅慧萍，李 杰

（1.上海交通大學 船舶海洋與建筑工程學院，上海200240；2.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240）

0 引 言

螺旋槳空泡是一種汽化空泡，即水因降壓到飽和蒸汽壓力導致汽化，水汽通過界面，進入氣核并使之膨脹。螺旋槳槳葉上出現(xiàn)的空泡有渦空泡、泡狀空泡、片空泡和云霧空泡4 種類型。螺旋槳空泡會對螺旋槳的水動力性能和噪聲性能產(chǎn)生重要影響。準確的數(shù)值模擬方法可以對螺旋槳性能進行預報，并為螺旋槳的設計提供理論依據(jù)。

近廿年來，國內(nèi)外對空化流提出的數(shù)值模擬方法有很多，大致可以分為勢流方法和Navier-Stokes方程求解及歐拉—拉格朗日混合方法三種。盡管作為一個水動力學熱點問題已經(jīng)被研究多年，但囿于目前的數(shù)值模擬準確性，仍然有眾多學者在空化模型及與空化相匹配的湍流模型，以及數(shù)值計算方法等方面開展相關的研究。尤其是對于螺旋槳梢渦及梢渦空泡的模擬，目前仍然是一個難點。德國波茲坦水池（SVA）曾于2011 年在漢堡主辦了第二屆船舶推進器國際研討會（SMP’11），該研討會采用盲算的方式進行比較分析，其中發(fā)布的第3 個空化案例（Case 2.3.3）在物理模型試驗中出現(xiàn)了顯而易見的梢渦空泡，數(shù)十家與會單位竟無一能算出。

國內(nèi)，韓寶玉等（2011）[1]針對梢渦空化問題，應用雷諾平均的N-S 方程并結合改進的VOF（volume fraction）模型，對橢圓水翼梢渦空化特性進行了研究，采用的湍流模型為代數(shù)雷諾應力模型。為了準確模擬梢渦渦核內(nèi)系統(tǒng)旋轉和流線曲率影響，將旋轉和曲率修正方法應用于湍流模型中。但是對水翼的水動力沒有進行定量的校核，只是給出了所謂消失空泡數(shù)。Shi 等（2015）[2]應用分離渦模型（DES）對螺旋槳尾流中梢渦流場分布進行了數(shù)值研究，但還沒有算過空化。Ji 等（2013）[3]采用大渦模擬（LES）和基于質(zhì)量輸運方程的空化模型預測了一個3 維扭曲水翼的非定?？栈?。辛公正 （2014）[4]采用數(shù)值RANS 方法模擬了水翼和螺旋槳的梢渦流動，并不求解空泡，而是尋找最小壓力系數(shù)與起始空泡數(shù)之間的關系。Liu 等（2016）[5]基于OpenFOAM 平臺，比較了線性與非線性湍流模型對螺旋槳梢渦空化的數(shù)值模擬能力，發(fā)現(xiàn)非線性k-ε 模型能夠更加準確地捕捉到梢渦空化現(xiàn)象。但沒有同時對螺旋槳的推力系數(shù)和扭矩系數(shù)進行校核。將螺旋槳空泡水動力性能計算做得比較好的有馮學梅等（2012）[6]，他們基于FLUENT 軟件，采用多塊結構化網(wǎng)格對業(yè)內(nèi)選作標槳的E779A 槳和PPTC 槳進行了均勻來流下的空泡數(shù)值模擬。

目前主流的空化計算方法是結合空化模型，求解兩相流的RANS 方程。歐拉—拉格朗日混合方法主要用于研究空化初生[7]。Hsiao 等（2005）[8]以水翼為研究對象，采用特定區(qū)域發(fā)放氣核的研究方法，考慮氣核尺寸分布進行梢渦空化初生噪聲尺度效應研究。真正與螺旋槳相關是在2008 年第27 屆海軍水動力學會議上，他們發(fā)表了敞水螺旋槳的梢渦空化初生尺度效應研究[9]。最近的研究成果是關于螺旋槳及氣體擴散作用對氣核分布的影響[10]。國內(nèi)熊鷹等（2013）[11]也采用這種方法對螺旋槳空化初生進行尺度效應研究，但均不直接模擬空泡。

本文基于螺旋槳空泡數(shù)值模擬的現(xiàn)狀，以PPTC 槳為研究對象，通過一種新穎的網(wǎng)格劃分方式，求解均勻來流下的螺旋槳空泡及其水動力。目標是同時求解出與試驗相符的空泡形狀，以及螺旋槳的推力和扭矩系數(shù)。

1 控制方程

采用多相流模型中的均質(zhì)混合流模型，此時氣液兩相流被當作包含微氣泡的單相連續(xù)介質(zhì)，通過引進氣相體積分數(shù)α 來定義混合相的密度ρm，僅求解一套控制方程?；诨旌厦芏鹊木|(zhì)混合流的連續(xù)方程為

動量方程為

式中：下標l、v 分別指液相和氣相；μm為混合流的動力粘性系數(shù)，定義與混合密度類似；μt,m為湍流引起的渦粘性系數(shù)，需要對混合流采用湍流模型進行求解。本文采用RNG k-ε 兩方程模型求解渦粘系數(shù)，該模型在求解旋轉流動和剪切流時有優(yōu)勢，壁面第一層網(wǎng)格節(jié)點的無量綱距離y+要求在30～100 之間。

關于蒸汽質(zhì)量分數(shù)f 的控制方程為

式中：Re和Rc分別為由汽化和凝結引起的相變率；f 按下式定義：

根據(jù)對（5）式右端源項的建模方式，可以區(qū)分出不同的空化模型。常用的有Zwart-Gerber-Belamri模型和Schnerr-Sauer 模型。本文采用Zwart 模型：

當p≤pv時，

當p＞pv時，

式中：pv為水的飽和蒸氣壓，與溫度相關；氣核半徑RB默認取值1 μm；αnuc為水中所含氣核的體積分數(shù)，默認取值5×10-4；Ce和Cc分別為汽化系數(shù)和凝結系數(shù)，默認取值50 和0.01。

2 研究對象與網(wǎng)格

PPTC 槳是指德國波茲坦水池（SVA）設計的可調(diào)螺距槳VP1304，也是第二屆和第四屆船舶推進器國際研討會（SMP’11 和SMP’15）的算例槳，其在均流和斜流情形下的敞水性能、伴流分布、空化性能及壓力脈動等試驗數(shù)據(jù)通過互聯(lián)網(wǎng)為全世界的水動力學研究人員所共享。本文即以SMP’11 的算例槳PPTC’11 為研究對象，對研討會上發(fā)布的空化案例Case 2.3 進行數(shù)值模擬，以校驗空化模型及相應的空化流場計算方法。

圖1 研究對象與計算域Fig.1 Computed object and computational field

PPTC 為直徑D=250 mm 的5 葉右旋可調(diào)螺距槳，這使得槳葉在導邊和隨邊處與槳轂之間存在0.3 mm 的間隙，數(shù)值模擬時予以忽略。SMP’11 的空泡測試在SVA 的空泡筒K15A 中完成，軸無傾角。分別在15、20、25 r/s 三個轉速下對敞水性能進行了測試，以探討雷諾數(shù)對敞水性能的影響；敞水試驗在高壓時進行，以避免空泡。圖1 為SVA 提供的PPTC 槳模與空泡筒試驗段幾何模型。為了更好地模擬均勻定常來流下槳的水動力性能，按SMP’11 的要求，計算時將試驗段變換成等截面的圓柱段，其它如槳模安裝位置和試驗段長度保持不變。

依據(jù)流動最大雷諾數(shù)及湍流模型的適用范圍估算了壁面第一層網(wǎng)格節(jié)點的距離，約為2.4D‰，據(jù)此在近壁區(qū)進行邊界層網(wǎng)格加密。尤其是對槳葉梢部的梢渦脫出區(qū)域進行了區(qū)域劃分及域內(nèi)網(wǎng)格加密，以保證渦核區(qū)域內(nèi)徑向15 個網(wǎng)格節(jié)點的密度（如圖2 所示）。生成全流道計算域網(wǎng)格，共8 560 925 個混合型單元。

圖2 梢部網(wǎng)格加密Fig.2 Local grid refinement of the tip vicinity

3 無空化流場計算

PPTC 槳被設計專門用于產(chǎn)生梢渦，并展開廣泛的模型試驗，用以校驗螺旋槳的各種勢流或粘流分析程序。在無空化情形下，SVA 在空泡筒對3 個轉速下的PPTC 敞水性能進行了測試。本章選擇3條敞水性能曲線中進速系數(shù)等于1 附近的3 個試驗點進行數(shù)值模擬。計算條件同試驗條件（見表1），表中pA為環(huán)境壓力。旨在獲得可靠的空化流場初始解，并對網(wǎng)格劃分及計算方法進行無空化計算校驗。

表1 無空化計算條件Tab.1 Computational conditions w/o CAV

計算在單相流框架下進行。湍流模型采用旋轉流適用的RNG k-ε 模型，近壁區(qū)的處理采用標準壁面函數(shù)。壓力速度耦合采用PISO（Pressure-Implicit with Splitting of Operators）算法，即隱式算子分割算法，該算法適用于高傾斜度網(wǎng)格?？臻g離散均采用二階以上高精度離散方式：梯度擴散項采用基于單元體的格林—高斯離散；壓力采用高旋流適用的PRESTO 格式；動量方程和湍流方程均采用適用于四面體網(wǎng)格的二階迎風格式。表2 給出了推力系數(shù)KT和扭矩系數(shù)KQ的計算結果。由表可見：KT和KQ的計算誤差在3.0%～5.4%以內(nèi)；隨著轉速提高（n=15～25 r/s），雷諾數(shù)增大（Re0.7=1.0×106～1.6×106，以0.7R 半徑處的合速度及槳葉弦長為特征量），計算精度略有降低。

表2 無空化計算結果與誤差分析Tab.2 Results and error analysis w/o CAV

為了準確定義并描述旋轉的流場，需要對漩渦結構進行定義。關于漩渦的判別法則有許多，比較流行的是渦量準則、Q 準則、Δ 準則和λ2準則。以上這些判別方式，都定義了一個可以逐點計算數(shù)值的函數(shù)，然后判定每一點位于漩渦內(nèi)或者漩渦外。盡管判定方式不少，但是對渦進行準確而嚴格的定義卻是較為困難的。通常，漩渦被認為是渦量較高的區(qū)域，但是對該區(qū)域的準確閾值卻沒有嚴格的定義，隨著流動性質(zhì)的不同，閾值的判定將會嚴重影響結果。而且，在沒有漩渦的強剪切流動中渦量也可能較大，所以，利用渦量準則來判定漩渦是不夠準確的。利用流線來定義渦并進行描述也是不可靠的，因為即使在簡單的勻速平移運動中，流線的拓撲結構也容易被改變。而后三種方法（Q 準則、Δ 準則和λ2準則）則在漩渦判斷的方式上具有等價性。本文采用Q 準則對漩渦進行定義和直觀的表達。

Q 準則由Hunt 等在1988 年提出[12]，它定義流場中速度梯度張量▽V 的第二不變量Q 具有正值的區(qū)域為漩渦。另外，它要求漩渦區(qū)域的壓強低于周圍的壓強。Q 定義為：

式中：Ω 為渦量，S 為變形率，皆為二階張量。符號“||||”表示張量的二范數(shù)。S 和Ω 分別代表了流場中一點的變形和旋轉。Q 準則反映了流場中一個流體微團旋轉和變形之間的一種平衡，Q＞0 則反映了旋轉在流動中占據(jù)統(tǒng)治地位，同時也體現(xiàn)了該旋轉區(qū)域的壓力為附近區(qū)域的最小值（極小值）。圖3 給出了以轉速為特征量的無量綱等值面圖。由圖可見，采取本文所述網(wǎng)格生成方法得到的梢渦流場較之采用傳統(tǒng)的不對梢渦進行區(qū)域劃分的網(wǎng)格計算結果要顯著得多，所捕獲的梢渦長度要長得多：梢渦從導邊后緣拖出，并向螺旋槳后方移動并逐漸耗散；當值增大，梢渦區(qū)域向渦核中心線收縮。

圖3 等值面Fig.3 ISO surfaces of

4 空化流場計算

螺旋槳空化性能預報對于分析螺旋槳的設計工況和非設計工況都意義重大。SMP’11 對工作狀態(tài)下的螺旋槳空化性能進行了觀測，均在轉速n=25 r/s 下進行，按不同的推力系數(shù)和空化數(shù)分為3 個案例，分別標注為Case 2.3.1、Case 2.3.2 和Case 2.3.3，試驗條件見表3。其中，Case 2.3.1 和Case 2.3.3 為非設計工況，Case 2.3.2 為設計工況。本章將對Case 2.3.1 和Case 2.3.2 進行數(shù)值模擬。計算條件同試驗條件。

表3 空化計算條件Tab.3 Computational conditions with CAV

表3 中σn為以旋轉線速度為特征速度的空化數(shù)，定義如下：由此可以換算出環(huán)境壓力pA。

先計算無空化情形，此時環(huán)境壓力可以取作一個大氣壓。校驗推力系數(shù)（扭矩系數(shù)缺失），并作為空化流場的初始解。然后將計算切換至均質(zhì)混合流框架下，改變環(huán)境壓力后，開啟空化模型并設置空化參數(shù)。此時，與無空化計算相比，增添了一個氣相體積分數(shù)的控制方程。該方程的空間離散選用QUICK（Quadratic upwind interpolation）格式。此格式用于六面體網(wǎng)格具有三階精度，用于四面體網(wǎng)格或混合網(wǎng)格時只具有二階精度。其他在無空化計算時已有的設置則保持不變。

4.1 算例1（Case 2.3.1）

該算例屬于非設計工況，載荷偏高情形。先進行無空化計算（此時pA=101.325 kPa），將計算得到的推力系數(shù)與試驗值進行比較。由表4 可見，推力系數(shù)的計算誤差為2.9%。重新設置環(huán)境壓力（此時pA=42.207 kPa），開啟空化計算。由表4 可見，空化將引起推力和扭矩系數(shù)下降；KT和KQ的計算誤差均在3.0%以內(nèi)。

表4 計算結果與誤差分析（Case 2.3.1）Tab.4 Results and error analysis for Case 2.3.1

空泡形態(tài)一般采用汽相體積分數(shù)琢的等值面來表示（0<琢≤1），因此空泡的大小長短與琢的取值相關。圖4 給出了空泡形態(tài)的計算結果（琢=0.2）與試驗結果（EFD）對比。其中，“SS”表示吸力面；“PS”表示壓力面；“SVSS”表示吸力面?zhèn)纫晥D。圖4（a）的試驗結果顯示，此工況下，除了發(fā)生顯著的梢渦空化外，吸力面葉根處也發(fā)生了片空化。圖4（b）～（d）的3 個視圖全方位地展示了梢渦空泡、吸力面葉根處的片空泡及轂渦空泡，計算與試驗符合較好。

圖4 空泡形態(tài)計算與試驗比較Fig.4 Observed and calculated cavitation(Case 2.3.1)

圖5 等值面Fig.5 ISO surfaces of (Case 2.3.1)

4.2 算例2（Case 2.3.2）

此工況為設計工況，中等載荷情形。先進行無空化計算（此時pA=101.325 kPa），將計算得到的推力系數(shù)與試驗值進行比較。由表5 可見，推力系數(shù)的計算誤差為4.5%。重新設置環(huán)境壓力（此時pA=30.528 kPa），開啟空化計算。由表5 可見：空化將引起推力和扭矩系數(shù)下降；KT的計算誤差為2.4%，KQ為4.5%。

表5 計算結果與誤差分析（Case 2.3.2）Tab.5 Results and error analysis for Case 2.3.2

圖6 給出了空泡形態(tài)的計算結果與試驗結果對比。圖6（a）的試驗結果顯示，此工況下的梢渦空泡在靠近隨邊處發(fā)生并拖出，吸力面葉根處發(fā)生了較大面積的片空化，并伴隨有泡狀空泡。圖6（b）～（d）分3 個視圖給出了數(shù)值模擬結果（琢=0.2）。由圖可見，除了梢渦空泡、吸力面葉根處的片空泡，壓力面葉根處也發(fā)生了較小的片空泡，計算與試驗符合較好。此工況下，沒有出現(xiàn)轂渦空泡。

圖6 空泡形態(tài)計算與試驗比較Fig.6 Observed and calculated cavitation(Case 2.3.2)

圖7 等值面Fig.7 ISO surfaces of (Case 2.3.2)

4.3 水中含氣率對空化計算的影響

數(shù)值實驗表明，當水中氣核體積分數(shù)琢nuc在一定范圍內(nèi)取不同值時，對空泡形態(tài)的影響并不顯著，卻會影響宏觀量—推力和扭矩系數(shù)的計算。表6 給出了算例Case 2.3.2 中琢nuc對KT和KQ的影響。由表可見，琢nuc越大，KT和KQ有降低的趨勢，從而更接近試驗值，但此時的計算容易發(fā)散。如何讓琢nuc的取值與模型試驗中的水中含氣量等環(huán)境變量相關聯(lián)，或者說，進一步探究或改良現(xiàn)有空化模型是下一步研究需要解決的問題。

表6 含氣率對KT 和KQ 的影響Tab.6 Effects of αnuc on KT&KQ

5 結 語

本文采用梢渦區(qū)域劃分的網(wǎng)格形式，對PPTC 槳進行了梢渦及梢渦空化模擬，在重載及中載工況的空化模擬中取得了初步的成效。通過數(shù)值模擬能得到顯著的梢渦及梢渦空泡，同時獲得較為準確的推力和扭矩系數(shù)。研究表明：梢渦空泡的捕捉很大程度上取決于梢渦的捕捉，梢渦空化只發(fā)生在渦量較大、壓力較低的渦核區(qū)域。