基于網格自適應技術的PPTC槳梢渦空泡數值模擬

陳 影，余 龍

（1.上海交通大學海洋工程國家重點實驗室；船舶海洋與建筑工程學院，上海 200240；2.高新船舶與深海開發裝備協同創新中心（船海協創中心），上海 200240）

0 引 言

螺旋槳空化問題一直是螺旋槳水動力研究中的熱點。目前研究螺旋槳空泡的手段主要有實驗測量和數值模擬兩種。由于空化模擬測試成本較高，大部分的研究都是基于數值模擬方法。螺旋槳的空泡主要分為片空泡、梢渦空泡和轂渦空泡等，目前許多國內外學者已經解決了槳葉表面片空泡的數值預測，但是關于梢渦空泡和轂渦空泡的研究尚不充分。

本文采用的螺旋槳為德國波茨坦水池（SVA）的模型槳VP1304，該槳作為基準螺旋槳于第二屆和第四屆國際船舶推進器專題研討會（smp’11和smp’15）被數十家單位計算，有大量的試驗數據和數值模擬結果[1-2]。Morgut[3]使用Zwart、FCM和Kunz的三個修正之后的傳輸模型對PPTC螺旋槳進行了CFD模擬，包括均勻流條件和傾斜軸條件；Gaggero[4]通過RANS方法基于OpenFOAM 開源平臺計算了PPTC螺旋槳在斜流中的無空化和空化非穩態性能；陳凱杰[5]基于OpenFOAM 開源平臺對PPTC 螺旋槳在全濕流和空化流條件下采用RANS和DES兩種方法進行了數值模擬。盡管這些文獻成功地驗證了PPTC螺旋槳的性能參數以及槳葉上片空泡的分布特性，但是無法模擬出梢渦空泡，尤其是脫出至尾流區域的梢渦空泡的模擬。

由于螺旋槳梢渦空泡尺度非常小，故梢渦空泡的數值模擬對網格分辨率要求很高，目前主流的梢渦空泡計算方法是根據初步計算確定梢渦空泡的大致范圍，然后確定控制體的螺距和半徑，接下來加密控制體內的網格單元。劉登成[6]采用方塊截面的控制體，胡健[7]和Yilmaz[8]采用的是螺旋管幾何體加密。雖然控制體加密的方法可以適當延長梢渦空泡，但是捕獲到的梢渦空泡均呈現出斷裂的現象，有一定的缺陷。故已有研究者試圖將自適應網格技術應用于梢渦空泡的模擬，以期得到更好的效果。

應用網格自適應技術的關鍵在于加密標準和加密區域單元尺寸的設置。加密標準直接影響到網格自適應的區域，應盡量保證加密區域與需要提高物理分辨率的區域重合，而在物理量變化較為平緩的區域保持不變或者適當地粗化網格，實現精準加密和提高計算精度，同時控制計算成本。Yilmaz[9]在模擬INSEAN E779AA 空化流時采用絕對壓強作為加密標準，當絕對壓力小于10 000 Pa（比飽和蒸氣壓高的壓力值）時，進行網格單元加密，取得了較好的效果，但網格增加了近4 倍。Eskilsson 等[10]將自適應網格技術應用于NACA 0015二維水翼的空化流模擬，并使用了4種不同的加密準則：速度u、渦量Q、體積分數α和壓強P，結果表明使用體積分數α最能精準地加密空化區域，而使用Q準則不僅能加密空泡區域，還加密了前緣和后緣區域。關于單元尺寸的控制，Kuiper（1981）[11]對V螺旋槳在J=0.3,0.4,0.5時梢渦空化的測量與研究中總結了空泡數和空泡半徑之間的經驗關系，發現空泡初生時，每個氣泡的最小半徑始終約為0.25 mm。故加密后的單元尺寸應該不高于0.25 mm，否則不能較好地捕獲到梢渦空泡。

本文采用了自適應網格技術研究梢渦空泡，使用Schnerr-Sauer空泡模型和SSTK-Ω湍流模型，采用體積分數α作為判據，提出一種高效率的網格劃分策略，實現空化區域的精準加密，選取了PPTC 槳作為計算對象，更精細地模擬了梢渦空泡。

1 數學模型與研究方法

1.1 控制方程

在數值模擬中，將空化流處理為包含液體和蒸氣的兩相流，采用修正之后的RANS 方程來求解。修正之后的連續性方程和動量方程如下：

式中：下標m表示多相流，ρm為多相流的密度，其余符號與常規RANS方程一致。

1.2 湍流模型

式中的經驗封閉系數和其余輔助方程可見文獻[12]，此處不再贅述。式（4）中的F1表示混合函數，其在近壁處邊界層上等于1，在遠離壁面處等于0，將K-Ω模型和K-ε模型結合在一起。

1.3 空泡模型

上述方程中的ρm和μm由氣相體積分數α確定，定義如下：

同樣地，需要補充一個關于α的輸運方程以使方程組封閉，即空泡模型，本文采用Schnerr-Sauer空泡模型[13]。Schnerr-Sauer 空泡模型在Rayleigh-Plesset 空泡模型基礎上，忽略了氣泡生長加速效應、粘性效應和表面張力效應。關于α的輸運方程為

1.4 自適應網格技術

自適應網格技術（Adaptive Mesh Refinement,AMR）是指基于計算結果根據自適應網格標準在計算過程中不斷細化網格單元，從而為某些物理值變化特別劇烈的區域提供足夠高的網格分辨率，在提高計算精度的同時又保證計算效率。對六面體網格來說，一個加密級別意味著1個父單元細化成8個子單元。

本文采用體積分數α作為控制標準，當α=0或1時，單元位于氣相內或者液相內，不進行加密；當α處于0和1之間時，相界面穿過該單元，進行加密。這樣就可以有效地自動加密空泡區域，在提高精度的同時又不會使網格數量過大。

2 研究對象與計算域設置

計算中采用SVA 提供的PPTC 螺旋槳幾何模型，該槳螺距可調，故在葉片和槳轂之間有很小的縫隙，模擬中予以忽略，其試驗模型和幾何參數分別見圖1 和表1。本文對研討會smp’11 上發布的空化案例Case2.3[14-15]進行數值模擬，以校驗無空化流場及空化流場的計算方法。

表1 PPTC模型槳的幾何參數Tab.1 Main particulars of PPTC propeller

圖1 PPTC槳試驗模型Fig.1 Model of PPTC propeller

關于螺旋槳的進速系數、推力系數、扭矩系數、敞水效率及空泡數的定義為

3 無空化流場計算與網格無關性分析

3.1 計算域設置

計算域（圖2）分為旋轉域和靜止域，采用速度入口和壓力出口條件，靜止域采用了smp’11[1]提供的空化流場試驗段，旋轉域直徑為1.2D。無空化流的計算采用多重參考系方法計算其穩態性能，空化計算采用滑移網格法計算非定常性能。

圖2 計算域及邊界條件設置（藍色區域為旋轉域）Fig.2 Computational domain and boundary condition settings for PPTC

3.2 網格劃分

本研究采用切割體網格劃分技術來劃分網格。在劃分網格時，整個計算域網格單元基準尺寸設為100 mm，對旋轉區域網格加密，加密區尺寸為基準值的8%。對葉片邊緣線網格以及葉片面網格進行局部加密，葉片邊緣附近最小網格尺寸為基準值的0.25%，槳葉表面最小網格尺寸為基準值的1%。螺旋槳葉片采用5層棱柱層，棱柱層總厚度為基準值的0.2%。網格（圖3）的數量為114萬。觀察計算完成后的Y+分布直方圖（圖3）可知，大部分網格單元的Y+小于1，基本上所有網格單元的Y+均小于5，這說明近壁單元基本上都位于粘性子層內。為了驗證網格無關性，本文設置了三種網格密度：64 萬、114萬和250萬。

圖3 114萬計算網格示意圖與Y+分布直方圖Fig.3 Computational mesh on the PPTC propeller and wall Y+distribution histogram

3.3 網格無關性分析

首先通過對三個進速系數（0.6928，1.2621，1.4944）下的無空化流場計算來驗證網格無關性，轉速均為15 r/s。表2 呈現了粗糙、中等、精細網格在三個進速系數下的KT和10KQ數值模擬結果以及網格增加時KT和10KQ的變化率。從該表中可以看出，變化率最高值為3.54%；而且隨著網格單元數量的增加，大部分計算結果呈現出降低的趨勢，只有J=0.6928 情況下KT的變化率有少許增加，但該變化率很小，未超過1%。總體而言，采用114 萬網格已經滿足計算結果對網格的無關性，故將此網格作為后文無空化流和空化流模擬的計算基準網格（G1）。

表2 推力系數和扭矩系數的網格無關性分析Tab.2 Grid independency analysis for the thrust coefficient and the torque coefficient

3.4 水動力性能參數分析

計算條件及計算結果如表3 所示，采用基準網格G1 進行其余工況的計算，表中ε表示計算結果和試驗結果之間的相對誤差。大部分工況下KT和10KQ的誤差小于1%，所有工況的誤差小于3%。圖4表明，計算值與試驗值吻合，說明了數值模擬的可靠性，可為進一步的空化模擬提供參考。

圖4 PPTC螺旋槳敞水特征曲線Fig.4 Propeller performance diagrams

表3 PPTC螺旋槳敞水性能計算結果及與試驗結果對比（CFD：數值模擬；EFD：模型試驗）Tab.3 Comparison of the open water propeller performance coefficients between measured and computed

4 空化流場模擬

4.1 計算參數設置

空化流場模擬時依賴初始條件的設置，故在進行空化流場模擬之前，先計算該空化流場對應的敞水情況，計算收斂后導出速度場和壓力場數據，然后將其導入至空化模擬中作為初始場，使計算快速收斂。將模型改為隱式非定常條件，時間步長設為一步旋轉1°的時間，時間離散格式為二階，計算參數與試驗[14]（見表4）保持一致，采用基準網格（G1）進行空化流場的模擬。

表4 PPTC槳空化流場模擬參數設置Tab.4 Experimental and computational conditions for cavitation flow

4.2 空化流場與敞水流場推力系數對比

表5 為三個案例中CFD 結果與EFD 結果之間推力系數的對比。整體而言，數值模擬的推力系數和試驗值較為接近，略低于試驗值；Case2.3.1 的無空化流和空化流的模擬效果最好，誤差較小；在Case2.3.2中，空化流的誤差偏高，達到5%左右；在Case2.3.3中，無空化流的誤差略大。

此外，從表5 還可以看出無空化流和空化流之間推力系數的差異，對每一個案例，空化的存在都會使螺旋槳的推力系數減小。特別是在Case2.3.2 和Case2.3.3 中，這一推力損失現象很突出，損失值達到了14%以上。

表5 無空化流和有空化流推力系數對比Tab.5 Comparison of the thrust coefficients between cavitation flow and no-cavitation flow

4.3 空泡分布情況

由圖5 可見，總體而言，三個案例的空泡分布區域與試驗較為接近，但細節上有所差別。在Case2.3.1 中，r＞0.95R范圍內的片空泡與試驗保持一致，接近槳轂的葉根區域空泡分布也與試驗很接近；但在導邊附近的0.5R至0.9R范圍內，數值模擬結果出現了試驗結果沒有的片空泡。在Case2.3.2中，數值模擬在導邊附近捕獲到部分空泡，與試驗不符；葉根區域的片空泡分布范圍略大于試驗。在Case2.3.3 中，空泡從吸力面轉移至壓力面上，數值模擬沒有捕捉到葉根區域的片空泡，且壓力面上導邊附近的空泡范圍略小于試驗。另外，三個案例的數值模擬均未捕捉到梢渦空泡，這是因為網格分辨率不足造成的。

圖5 空泡分布形態對比（數值模擬中用α=0.2等值面表征空泡）Fig.5 Comparison between the computed cavitation and the experimental cavitation

5 梢渦空泡

5.1 網格自適應加密

為了捕獲尾流中的梢渦空泡，需提高該處的網格分辨率。本文采用三套網格（圖6）來模擬梢渦空泡。第一套網格就是前文所述的基準網格（G1），網格單元數量為114萬。第二套網格（G2）在G1基礎上應用螺旋管幾何加密，螺旋管幾何如圖6 所示，單元數量為415 萬；其中螺旋管的直徑為10 mm，其螺距和半徑減少量從G1 計算結果中獲得，加密基本尺寸為0.5 mm。第三套網格（G3）是以G2 作為初始網格在計算中應用本文采用的網格自適應技術后生成的最終網格，數量為638萬；當氣體體積分數α在0至1之間時自動加密單元，加密級別取1，加密后單元基本尺寸變為0.25 mm，網格優化貫穿整個計算過程，每5個時間步應用一次網格自適應。若不應用網格自適應方法，直接將螺旋管內單元尺寸設為0.25 mm，則網格量高達2192 萬，故本文中的網格劃分策略可顯著減少網格單元數量，提高計算效率。文中計算梢渦空泡的具體流程圖如圖7所示。

圖6 螺旋管及計算網格Fig.6 Spiral tube geometry and three kinds of computational mesh

圖7 梢渦空泡計算流程圖Fig.7 Flow chart of tip vortex cavitation simulation

5.2 推力系數對比

表6 為三種網格劃分策略的推力系數計算結果的對比，其中KTG1表示采用網格G1 計算得到的推力系數，KTG2、KTG3同理。三套網格計算的結果僅有0.6%以內的微小差異，幾乎可忽略。這三套網格之間最大的不同在于對梢渦脫出區域的處理不同，對推力系數的影響很小，在正常的波動范圍之內。

表6 三種網格劃分策略推力系數對比Tab.6 Comparison between CFD and EFD for three mesh methods

5.3 空泡分布形態對比

圖8展示了三種網格劃分策略的空泡分布形態，槳葉上的片空泡基本沒有變化，但采用螺旋管與網格自適應技術結合的網格劃分策略可以明顯改善梢渦空泡和轂渦空泡的模擬，特別是尾流區域中的梢渦空泡延長效果顯著。

圖8 梢渦空泡的延長效果對比圖Fig.8 Improvement of tip vortex cavitation extension

5.4 梢渦空泡細節對比

由圖9（a）可見，梢渦空泡呈現出卷起的現象，這一現象是由于渦流強度降低或壓力增加，導致片空泡和梢渦空泡產生相互作用而引起的，卷起時產生一定規律的節點，節點處的空化渦管直徑減小。僅采用螺旋管加密的G2網格模擬出的梢渦（圖8）在節點處由于網格分辨率不足，呈現出斷裂的現象；但是觀察圖9（b）可知，應用網格自適應技術后成功地模擬了渦管卷起現象，更精確地獲得了流場的空化模式。

圖9 Case2.3.1梢渦空泡放大對比圖Fig.9 Comparison of the tip vortex cavitation roll-up phenomena between EFD and CFD

5.5 漩渦結構圖

采用Q準則來表達漩渦，當Q＞0時表示旋轉在流動中占主要地位，其值越大表示渦的強度越大。圖10第二列展示了G3計算結果，用Q為10 000的等值面來表示漩渦。與G1計算結果（圖10第一列）相比，捕獲的梢渦長度要長得多，而且呈現出更為光滑的管狀結構，渦管的直徑也更加均勻，但是到后面未加密區域時，渦管直徑變大，并且很快耗散，若想延長捕獲的渦管，可通過增加螺旋管的長度來實現。

圖10 梢渦與梢渦空泡對比圖Fig.10 Comparison between tip vortex and tip vortex cavitation

觀察Q=10 000的等值面與α=0.2的梢渦空泡等值面（圖10第三列）可以發現，兩者形態十分接近，但是可以觀察到Q的等值面范圍遠遠大于α的范圍，特別是槳葉表面上，此處可證明采用Q準則作為網格自適應的判斷標準時會進行許多不必要的單元加密。若僅想探究螺旋槳的梢渦空泡，從加密精準度來看，氣體體積分數α是一個更好的選擇。

圖11展示了兩個案例x=0.3R截面處的渦量分布圖，左半邊為采用G1計算的結果，右半邊為采用G3 計算的結果。可以看到，兩種網格的渦量分布圖均為周期性分布，周期為5，對應于螺旋槳葉數。但是采用G3網格計算可以獲得更精細的渦量分布，由于G1網格在梢渦處分辨率不夠高，所以捕捉不到渦核處的最大渦量，而應用螺旋管精細加密和網格自適應的梢渦則顯示出更大的渦強。此外，結合圖12 可以發現渦量最大值并非出現在梢渦的中心，在向梢渦中心靠近時，渦量強度呈現出先增加后減少的變化趨勢。這是由于空泡的作用，當空化流中的壓力低于飽和蒸氣壓時，就會產生空泡，使壓力不會再降低，因此在空泡內部壓力相對比較均勻，故渦量最大值并非出現在中心。

圖11 渦量分布云圖對比（x=0.3R處）Fig.11 Vorticity comparison between different meshes

圖12 Q準則分布圖（Case2.3.1：x=0.3R，r=0.96R處）Fig.12 Q criterion distribution across the tip vortex(Case2.3.1:x=0.3R,r=0.96R)

5.6 速度分布

由圖13可知G3網格獲得了更精細的速度分布。G3結果在下圖紅色方框中呈現出更大的速度梯度，在遠離槳轂側的低速區和紅色的高速區的渦相互一一對應，并且耗散得較慢。

圖13 速度分布云圖Fig.13 Velocity distribution

6 結 論

本文通過RANS 方法對PPTC 螺旋槳無空化和有空化流場的模擬，提出一種新的網格劃分策略，并將其應用至梢渦空泡的模擬中，得到以下結論：

（1）通過RANS 方法采用未加密的網格計算無空化流場和空化流場，水動力系數誤差最高為6.26%，已達到較高的精度；模擬得到的槳葉片空泡分布區域與試驗基本相符，但是未能捕獲梢渦空泡。

（2）本文提出的采用螺旋管加密與基于氣體體積分數的網格自適應方法相結合的方法，可實現空化區域的精準加密，更高效地模擬出梢渦空泡。該方法有望廣泛應用于其它螺旋槳梢渦以及梢渦空泡的數值模擬中，可以用較少的網格預報槳葉片空泡以及梢渦空泡，在保證精度的同時控制計算成本。

（3）成功將新的網格劃分策略應用于PPTC螺旋槳的空化模擬中，延長了捕獲的梢渦空泡長度，還更好地呈現了梢渦空泡上卷的細節特征。之后也進行了尾流場中梢渦和速度分布相關分析，捕獲到顯著的梢渦，同時發現空化會使渦量最大區域從中心轉移。