劉芳遠(yuǎn)，傅慧萍，李 杰

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

0 引 言

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

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

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

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

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

1 控制方程

采用多相流模型中的均質(zhì)混合流模型，此時(shí)氣液兩相流被當(dāng)作包含微氣泡的單相連續(xù)介質(zhì)，通過引進(jìn)氣相體積分?jǐn)?shù)α 來定義混合相的密度ρm，僅求解一套控制方程。基于混合密度的均質(zhì)混合流的連續(xù)方程為

動(dòng)量方程為

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

關(guān)于蒸汽質(zhì)量分?jǐn)?shù)f 的控制方程為

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

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

當(dāng)p≤pv時(shí)，

當(dāng)p＞pv時(shí)，

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

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

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

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

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

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

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

3 無空化流場計(jì)算

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

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

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

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

為了準(zhǔn)確定義并描述旋轉(zhuǎn)的流場，需要對漩渦結(jié)構(gòu)進(jìn)行定義。關(guān)于漩渦的判別法則有許多，比較流行的是渦量準(zhǔn)則、Q 準(zhǔn)則、Δ 準(zhǔn)則和λ2準(zhǔn)則。以上這些判別方式，都定義了一個(gè)可以逐點(diǎn)計(jì)算數(shù)值的函數(shù)，然后判定每一點(diǎn)位于漩渦內(nèi)或者漩渦外。盡管判定方式不少，但是對渦進(jìn)行準(zhǔn)確而嚴(yán)格的定義卻是較為困難的。通常，漩渦被認(rèn)為是渦量較高的區(qū)域，但是對該區(qū)域的準(zhǔn)確閾值卻沒有嚴(yán)格的定義，隨著流動(dòng)性質(zhì)的不同，閾值的判定將會(huì)嚴(yán)重影響結(jié)果。而且，在沒有漩渦的強(qiáng)剪切流動(dòng)中渦量也可能較大，所以，利用渦量準(zhǔn)則來判定漩渦是不夠準(zhǔn)確的。利用流線來定義渦并進(jìn)行描述也是不可靠的，因?yàn)榧词乖诤唵蔚膭蛩倨揭七\(yùn)動(dòng)中，流線的拓?fù)浣Y(jié)構(gòu)也容易被改變。而后三種方法（Q 準(zhǔn)則、Δ 準(zhǔn)則和λ2準(zhǔn)則）則在漩渦判斷的方式上具有等價(jià)性。本文采用Q 準(zhǔn)則對漩渦進(jìn)行定義和直觀的表達(dá)。

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

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

圖3 等值面Fig.3 ISO surfaces of

4 空化流場計(jì)算

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

表3 空化計(jì)算條件Tab.3 Computational conditions with CAV

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

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

4.1 算例1（Case 2.3.1）

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

表4 計(jì)算結(jié)果與誤差分析（Case 2.3.1）Tab.4 Results and error analysis for Case 2.3.1

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

圖4 空泡形態(tài)計(jì)算與試驗(yàn)比較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）

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

表5 計(jì)算結(jié)果與誤差分析（Case 2.3.2）Tab.5 Results and error analysis for Case 2.3.2

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

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

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

4.3 水中含氣率對空化計(jì)算的影響

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

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

5 結(jié) 語

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