定常空化下船用螺旋槳推力損失機理數值分析

劉 謙 李子如 何朋朋 賀 偉*

(武漢理工大學高性能艦船技術教育部重點實驗室1) 武漢 430063) (武漢理工大學船海與能源動力工程學院2) 武漢 430063)

0 引 言

空化是一種包含相變、可壓縮性等特點的復雜流動現象.根據空化的發生位置及形態，螺旋槳空化一般可分為泡空化、片空化、超空化、渦空化和云空化五類.螺旋槳空化對螺旋槳的水動力性能及噪聲性能有深遠影響，而準確的數值模擬方法是探明這種影響的基礎，也可以為螺旋槳的設計提供理論依據.

早期針對螺旋槳空化問題的研究，主要采用基于勢流理論的邊界元方法預報螺旋槳的片狀空泡，后續研究開始采用基于黏流的計算流體力學(computational fluid dynamics，CFD)方法.韓寶玉等[1]應用雷諾平均的N-S方程并結合改進的VOF(volume of fluid)模型，研究了橢圓水翼梢渦的空化特性.劉志輝等[2]基于OpenFOAM平臺，比較了線性與非線性湍流模型對螺旋槳梢渦空化的數值模擬能力，發現非線性k-ε模型能更準確的模擬梢渦空化，但二者均未對研究對象的水動力性能進行定量校核.馮玉梅等[3]基于FLUENT軟件，采用多塊結構化網格對E779A槳和PPTC槳進行了均勻來流下的空泡數值模擬，其所預報的空泡形態、螺旋槳推力系數及扭矩系數均與試驗結果一致.劉登成等[4]采用Schnerr-Sauer空化模型，分析了梢渦與梢渦空泡的流動特征，指出空泡流中梢渦空泡區域的渦量的周向分布呈現雙峰特性.吳家鳴等[5]研究了空化條件下導管螺旋槳周圍流場及推力特性的變化，指出空泡因素就螺旋槳葉水動力特性而言，對吸力面及壓力面均有不可忽略的影響.劉芳遠等[6]采用Zwart-Gerber-Belamri空化模型，通過梢渦區域的劃分及網格加密對空化流場進行了數值模擬，所得螺旋槳推力系數和轉矩系數與試驗值相吻合.胡健等[7]采用大渦模擬(large eddy simulation，LES)方法和Schnerr-Sauer空化模型較好地預測了E779A螺旋槳的梢渦空化.胡洋等[8]基于RANS方法預測了斜流對槳葉載荷的影響.

針對螺旋槳空化問題的大多數研究集中在空泡形態預報上，而忽略了無空化流場與空化流場中螺旋槳推力的比較.文中以PPTC槳為例，通過比較無空化與空化流場中槳葉表面的壓力分布差異，探究空化造成螺旋槳推力下降的機理.

1 數學模型

采用多相流模型中的均質混合流模型，將汽液相視為同一介質，通過引進汽相體積分數α來定義混合相的密度ρm，使流體可通過一套方程求解.基于混合密度的均質混合流的連續方程、動量方程分別為

(1)

(2)

式中：ui為流體i方向速度；Fi為體積力；τij為黏性應力張量.介質為牛頓流體時定義為

(3)

式中：μ為動力黏性系數；δij為克羅內克爾符號(當i=j時δij等于1，否則其等于0).

式(2)中ρm為混合相的混合密度，為

ρm=αρv+(1-α)ρl

(4)

式中：ρv為汽相密度；ρl為液相密度；α為汽相體積分數.

蒸汽質量分數f的控制方程為

(5)

f=ρvα/ρm

(6)

選擇Schnerr-Sauer空化模型模擬空化流場，該模型采用氣泡數密度及氣泡半徑定義氣相體積分數，并且不考慮非冷凝氣體及湍流脈動對空化流的影響.

(7)

式中：pV為汽化壓力；p為當地壓力；RB為氣泡半徑，表示為

(8)

式中：e為氣泡數密度，此處選擇1013.

2 計算模型

選取VP1304螺旋槳(以下簡稱PPTC槳)為研究對象，其幾何圖形見圖1，相關參數見表1.

圖1 PPTC槳示意圖

表1 PPTC槳幾何參數

計算基于Fluent軟件在均質多相流的框架下進行.計算域參考空泡水筒試驗段的形狀及尺度設為10.3D×2.4D×2.4D的長方體，入口、出口到槳盤面距離分別為2.3D和8D.計算域入口設置為速度入口，出口設置為壓力出口；螺旋槳表面及計算域四周均設置為無滑移壁面條件；旋轉域采用MRF方法來模擬螺旋槳的旋轉運動，靜止域和旋轉域通過交界面來進行數據傳遞.計算域相關設置見圖2.

圖2 螺旋槳計算域

湍流模型采用旋轉流適用的RNGk-ε模型，空化模型如前所述選為Schnerr-Sauer模型，近壁面處理采用Non-Equilibrium Wall Functions.壓力速度耦合采用Coupled算法，梯度擴散項采用基于單元體的格林-高斯離散，壓力插值采用高旋流適用的PRESTO(pressure staggering option)；動量方程、氣相體積分數的控制方程和湍流方程均采用QUICK(quadratic upwind interpolation)格式.

采用混合網格劃分策略，靜止域采用六面體結構化網格，旋轉域采用四面體網格，在槳葉壁面附近進行網格加密.

表2 各方案的網格數量

3 結果與討論

3.1 網格不確定度分析

以進速系數J=1.019 3，空泡數σn=2.024為例進行網格不確定度分析，各參數定義如下.

進速系數

(9)

推力系數

(10)

扭矩系數

(11)

效率

(12)

空泡數

(13)

式中：Va為進速；螺旋槳轉速n設為24.987 r/s；流體密度ρ=997.44 kg/m3；T為螺旋槳總推力；Q為螺旋槳轉矩；空泡數的變化通過調整出口壓力p∞來實現，空化壓力取為pV=2 873 Pa.

表3為各網格下計算結果及對應的試驗結果，不確定度分析相關詳細計算過程略去不表，最終結果見表4～5.表中M表示對應的試驗結果由表4～5可知：推力系數kt及扭矩系數kq的計算結果均呈單調收斂，且|E|

表3 各方案數值結果及試驗結果

表4 推力系數kt驗證及確認

表5 轉矩系數kq驗證及確認

圖3 槳葉表面及流場局部網格示意

3.2 不同空泡數下空化流場分析

選取進速系數J=1.019 3的五個空泡數工況(σn=2.024、3.026、4.028、4.529、5.026)進行數值模擬，并與SVA實驗室提供的試驗結果進行比較，以進一步驗證前述網格劃分及數值方法的可行性.

表6為不同空泡數下的計算結果.由表6可知：kt和kq的計算誤差均在4%以內；其中，σn=3.026～5.026工況下的計算誤差均在3.2%左右，高于σn=2.024時0.1%左右的計算誤差.

表6 數值計算結果與試驗結果相對誤差比較

在數值計算中，空泡形態一般采用汽相體積分數α的等值面來表示(0<α≤1)，因此空泡的大小長短與α的取值相關.本文展示空泡形態時統一取α=0.2.圖4為σn=2.024工況下數值計算與試驗觀測得到的空泡形態對比.由圖4可知.

1) 數值計算可以在隨邊靠近葉梢處捕捉到較為明顯的梢渦空泡，但由于平均意義下RANS方法所固有的局限性以及當前網格密度稀疏導致數值計算無法捕捉到較大范圍的槳后梢渦；

2) 數值計算在槳葉根部可以捕捉到片狀空泡現象，但無法捕捉到如試驗結果所示片空泡破碎形成的泡狀空泡.其原因可能在于本文采用均質混合流模型和Schnerr-Sauer空化模型，將空泡過程中的蒸發、冷凝現象直接與流場壓力相關聯，對單個泡狀空泡的捕捉略顯不足.

圖4 螺旋槳在J=1.019 3,σn=2.024工況下的空泡形態

梢渦渦核處壓力低是梢渦空化形成的主要成因.本節還采用最為常用的Q準則對螺旋槳梢渦進行進一步的判斷分析.Q準則定義為

(14)

式中：‖Ω‖和‖S‖為渦量張量和應變率張量的二范數.Q準則反映了流場中流體微團旋轉和變形之間的平衡，當Q準則大于0時意味著旋轉在流動中占據統治地位.進一步根據轉速n對Q準則進行無量綱化，定義為

(15)

圖5 螺旋槳在J=1.019 3,σn=2.024工況下的空泡形態數值計算結果

綜上所述，雖然本文所用網格在捕捉槳后梢渦方面稍有遜色，但梢渦空泡對槳葉壓力分布的影響不大，故可以認為本文所用網格及計算方法就預報槳葉壓力分布及螺旋槳推力而言具有可行性.

3.3 不同進速系數下空化流場分析

保持σn=2.024不變，分別選取5個進速系數工況(J=0.573 9、0.787 5、1.019 3、1.268 0、1.408 3)進行數值計算.為便于比較，還對上述工況無空化情況進行了數值計算.表7為各工況下無空化及空化條件下的kt對比，相應的kt隨J的變化曲線見圖6.由圖6可知：發生空化時螺旋槳推力將出現下降；隨著進速系數增大，因空化發生導致螺旋槳推力下降的程度逐漸減小.

表7 無空化和空化條件下kt比較

圖6 無空化及空化條件下kt隨J的變化曲線

圖7為不同進速系數下無空化情況槳葉壓力系數分布情況.壓力系數基于旋轉線速度進行無量綱化，其定義為

(16)

式中：各物理量含義與式(7)、式(13)相同.

由圖7可知：J=0.573 9時，葉背負壓區在弦向具有一定的寬度，且從葉根順著導邊一直延伸到葉梢區域；隨著進速系數增大，葉背負壓谷值區域向導邊收縮，且葉片徑向中間區域的低壓區幅值及面積均呈減小趨勢；至J=1.408 3時，在葉面導邊開始出現負壓區，其谷值與低進速系數時葉背導邊的負壓谷值相當.該圖清晰的反應了隨著進速系數的增大，槳葉表面低壓谷值區域逐漸由吸力面向壓力面、葉梢向葉根的移動過程.

圖7 無空化流場下壓力面(左)和吸力面(右)壓力系數分布

圖8為以上三個進速系數工況對應的有空化情況下的槳葉表面壓力系數分布.圖9為對應的螺旋槳空泡形態數值計算結果，圖中黑色弧線標識了葉剖面徑向位置，其中J=0.573 9、1.019 3時，槳葉壓力面無空泡，故僅給出了吸力面一側的空泡形態.

圖8 空化流場下壓力面(左)和吸力面(右)壓力系數分布

當不考慮空化時，槳葉表面的低壓可以下降到汽化壓力以下，而當考慮空化時，空泡覆蓋區域的壓力最低只能下降至汽化壓力.其次，除了發生空化區域的壓力較無空化情況有所不同以外，未發生空化區域的壓力分布一定程度上也會受到空化的影響而發生改變.本節結合不同進速系數下的空泡形態，通過提取不同徑向位置處葉剖面的壓力系數分布并與無空泡情況比較，還進一步分析了不同進速系數工況下考慮空化時螺旋槳推力下降程度與以上兩點原因的關系.

圖9 各進速系數下螺旋槳空泡形態

由圖9a)可知：當J=0.573 9時，吸力面空泡覆蓋面積較大，從葉根順著導邊一直延伸到葉梢，包括整個外半徑區域(r/R≥0.8).圖10為J=0.573 9時有無空化情況下各葉剖面壓力系數分布，由圖10a)～10b)可知：與無空化情況相比，空化情況下吸力面在空化發生區域內始終維持汽化壓力(Cp=-σn)，由此導致槳葉在外半徑區域葉面與葉背兩側的壓力差更小，相應地推力也就會出現較明顯的下降；由圖10c)～10d)可知：在槳葉根部發生空化區域附近槳葉表面仍維持汽化壓力，而在未發生空化的壓力面，空化流場中弦向大部分區域的壓力均略低于無空化流場，靠近隨邊處的變化則相對復雜，且明顯受到吸力面空泡發生的影響.整體上看，考慮空化時槳葉葉根區域的兩側壓力差較無空化情況仍有所減小，但其下降幅度相較于葉梢區域明顯緩和.

圖10 J=0.573 9時有無空化情況下各葉剖面壓力系數分布

由圖9b)可知：當J=1.019 3時，吸力面空泡覆蓋面積較J=0.573 9工況大幅減小，僅在中等半徑處的導邊，以及葉根和葉梢靠近隨邊處有少量空泡.圖11為該進速系數工況有無空化情況下典型剖面壓力分布對比也大體與之對應.由圖11可知：

1) 在槳葉葉梢處，近隨邊發生空化區域(0.65≤x/c≤1.0)槳葉吸力面與壓力面的壓力較無空化情況均有一定的變化，且吸力面一側的變化明顯更大.

2) 在槳葉中等半徑處，始于導邊的小體積空泡對該位置剖面的壓力分布影響較小，僅空泡發生區域的壓力維持為汽化壓力，其他區域壓力與無空化情況基本相當.

3) 葉根僅在r/R<0.4的區域出現空泡.在r/R=0.35處，空泡區域(0.2≤x/c≤0.6)維持壓力造成的壓差減少量已明顯下降，在近隨邊區域(0.65≤x/c≤1.0)，壓力面壓力升高，吸力面壓力下降，此時壓力分布改變造成的壓差減小與空泡區域壓力無法進一步下降造成的影響處在同一量級.r/R=0.4處剖面壓力分布總體上也體現了葉根空泡發生的影響，但其變化造成的壓差減小量明顯要低于r/R=0.35剖面.

圖11 J=1.019 3時有無空化情況下各葉剖面壓力系數分布

圖12為進速系數J=1.408 3工況有無空化情況下典型剖面壓力分布對比.由圖12可知：在外半徑區域(r/R>0.75)已無空化出現，該區域槳葉兩側的壓力分布及壓力差與無空化情況相比無明顯變化；在中等半徑處(0.4

在葉根區域(r/R≤0.4)，無空化情況下，導邊側(0

圖12 J=1.408 3有無空化情況下各葉剖面壓力系數分布

與無空化情況相比，發生空化時其附近區域的壓力都會出現變化進而導致槳葉推力出現下降.在低進速系數工況下，發生空化區域壓力只能降低至汽化壓力導致壓差下降占主導；隨著進速系數的增大，空泡面積減小，空泡發生對附近區域壓力分布變化導致的壓差減小量所占比重逐漸增大.

4 結 論

1) 空化現象造成螺旋槳推力下降的兩個原因：①空化導致該區域壓力維持在汽化壓力左右，這相較于無空化流場中，阻礙了該區域負壓谷值的進一步降低；②空化導致未發生空化區域的壓力分布改變.

2) 當進速系數較低時(J=0.573 9)，空泡形態以覆蓋整個吸力面外半徑區域的片空泡為主，此時造成螺旋槳推力下降以上述第一個原因為主.而當進速系數逐漸增大(J=1.019 3、1.408 3)，槳葉空泡覆蓋面積迅速減少，空泡在槳葉表面即發生潰滅，空化區域維持汽化壓力對槳葉壓差的影響逐漸減小，造成螺旋槳推力下降的主要原因從上述第一點轉變為第二點.