非均勻流場螺旋槳空泡數(shù)值模擬

劉 恒，伍 銳,3，孫 碩

(1.上海船舶運輸科學研究所 航運技術(shù)與安全國家重點實驗室，上海 200135；2.上海船舶運輸科學研究所 航運技術(shù)交通行業(yè)重點實驗室，上海 200135；3.浙江大學 能源工程學院化工機械研究所，杭州 310027)

19世紀末，英國學者就開始注意到螺旋槳空化現(xiàn)象，空泡的發(fā)生導致螺旋槳工作于水汽中，螺旋槳推力下降，航速難以達到預期[1].空化其實是一種熱力學現(xiàn)象，環(huán)境壓力不變當溫度升高時，液體吸收熱量發(fā)生汽化；溫度下降時，水汽冷凝成液體狀態(tài).對于螺旋槳空化而言,溫度一定，槳葉高速旋轉(zhuǎn)，槳葉表面壓力低至飽和蒸汽壓力時，液體水分子汽化，水汽通過界面進入氣核使之膨脹發(fā)展為肉眼可見空泡；當槳葉表面壓力升高，空泡發(fā)生潰滅，釋放熱量，產(chǎn)生劇烈壓力脈沖.常見的螺旋槳空泡形態(tài)有片空泡、梢渦空泡、泡空泡、葉根空泡、云霧狀空泡及轂渦空泡等，其中梢渦空化起始最早，會導致螺旋槳輻射噪聲驟增，水下航行器臨界航速的高低完全取決于此.因此，如何延緩梢渦空泡的產(chǎn)生對于槳的設(shè)計者而言至關(guān)重要，而準確預測梢渦空化則是低噪聲螺旋槳設(shè)計的前提和基礎(chǔ).

相對于試驗而言，數(shù)值模擬能夠以較低的成本、更高的效率及能提供更加豐富的流場信息而成為研究螺旋槳空化流動的重要手段.空化數(shù)值模型大體分為三類：界面追蹤法[2]、全流場求解Navier-Stokes(N-S)方程均相流模型[3]及歐拉-拉格朗日混合方法[4].其中N-S方程全流場求解應用最為廣泛，該方法將液相和氣相均視為連續(xù)流體介質(zhì).朱志峰等[5]以E779a螺旋槳模型為研究對象，采用Singhal全空化空化模型，在勻流場中比較了標準k-ω模型、修正Renormalization-group (RNG)k-ε模型和大渦模擬方法(LES)大渦模型對空化后螺旋槳推力系數(shù)和蒸汽體積分數(shù)的影響.其計算結(jié)果表明：標準k-ω模型能適應質(zhì)量較低的非結(jié)構(gòu)網(wǎng)格，穩(wěn)定性較好，但對速度和湍流動能的預報精度不高；修正RNGk-ε模型對葉面隨邊處邊界層周圍黏性尾流產(chǎn)生的速度凹陷預報較為準確，對梢渦附近速度的預報也具有一定的精度.相比而言，LES在預報速度和湍流動能時比求解時均N-S方程的前兩種湍流模型更精確，但計算量更大，耗時更久.Liu等[6]同樣以E799a為研究對象，采用Schnerr-Sauer 空化模型，在勻流場中比較了標準k-ε湍流模型和非線性k-ε湍流模型對梢渦空泡的捕捉能力，其研究表明非線性k-ε湍流模型能夠有效預測梢渦空化.胡健等[7]采用LES和Schnerr-Sauer 空化模型，依然以E799a為研究對象，通過對螺旋槳葉梢附近處網(wǎng)格加密(網(wǎng)格單元數(shù)量近1 900萬)，成功捕捉到勻流場中螺旋槳梢渦空化，數(shù)值結(jié)果與實驗結(jié)果吻合良好.劉芳遠等[8]采用RNGk-ε湍流模型和Zwart-Gerber-Belamri(Z-G-B)空化模型，對勻流場中PPTC(Potsdam Propeller Test Case)槳進行了空化數(shù)值計算，有效地實現(xiàn)了梢渦捕捉，同時指出水中含氣率對推力系數(shù)和轉(zhuǎn)矩的影響比空泡形態(tài)更強.

螺旋槳一般工作于非均勻流場中，上述的研究均是在均勻流場中進行的，并未對更加復雜的非均勻流場進行分析.Ji等[9]在非均勻流條件下，采用SSTk-ω湍流模型和Z-G-B空化模型對E799a螺旋槳空泡性能進行了數(shù)值分析，研究表明其數(shù)值計算葉背片空泡形態(tài)、空化體積與試驗結(jié)果吻合良好，但并未捕捉到梢渦空化.傅慧萍等[10]以PPTC槳為研究對象，結(jié)合RNGk-ε湍流模型和Z-G-B空化模型探究了斜流中螺旋槳的空化形態(tài)，但是其片空泡、梢渦空泡計算結(jié)果與試驗值略有差異.

本文以某油船槳為研究對象，基于RANS(Reynolds-averaged Navier-Stokes)方程，采用Schnerr-Sauer 空化模型和可實現(xiàn)的k-ε兩層湍流模型，采用商業(yè)軟件STAR-CCM+求解非均勻流場中螺旋槳的空泡形態(tài)，將計算得到的各相位角下葉背片空泡正投影面積、葉背片空泡形態(tài)及梢渦空泡形態(tài)與試驗結(jié)果進行對比.目標是求解出各相位角下與試驗相符的空泡形態(tài)，建立快速、高效、準確的空泡數(shù)值預報方法.

1 控制方程和空泡模型

在氣液兩相流混合模型中，流場被認為由可變密度的單一流質(zhì)組成，忽略相間的速度滑移，氣相與液相具有相同的速度和壓力.在慣性直角坐標系下，定義坐標軸x為螺旋槳的旋轉(zhuǎn)軸，沿水流方向為正，見圖1.因此基于混合密度的均質(zhì)混合流的連續(xù)性和動量控制方程為

(1)

(2)

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

(3)

μm=μvαv+μlαl

(4)

式中：下標i和j分別為坐標方向；下標 l、v分別指液相和氣相；t為時間；ρm為混合流密度；u和p分別為速度和壓力；μm為混合相的動力黏性系數(shù)；μt,m為湍流黏性系數(shù)；α為相體積分數(shù).ρm和μm均通過體積加權(quán)平均后得到.采用可實現(xiàn)的k-ε湍流模型求解渦黏性系數(shù)，該模型適用于強旋轉(zhuǎn)流動和剪切流，文獻[2-5,7]采用該湍流模型得到了較為滿意的結(jié)果.

由于控制方程中增加體積分數(shù)這一未知量，因此需通過添加額外的方程進行封閉.相變過程可認為起始于溶解于水中的氣核，單位控制體積內(nèi)混合流的質(zhì)量始終保持不變，根據(jù)氣核的生長和潰滅速率可計算兩相間的質(zhì)量輸運率和體積分數(shù).基于Rayleigh-Plesset方程的Schnerr-Sauer空化模型能夠?qū)庖簝上嚅g的質(zhì)量輸運進行較好的描述[11].體積分數(shù)的控制方程為

(5)

(6)

(7)

式中：p0為飽和蒸氣壓，取值 3 850 Pa，試驗時水溫為28.4 ℃.隨著空化逐漸發(fā)展，單位體積內(nèi)汽泡個數(shù)N0=αln0，保證液相中汽泡密度始終保持不變.

當p

(8)

當p>p0時，氣體凝結(jié)，質(zhì)量輸運速率為

(9)

2 研究對象、網(wǎng)格劃分及可信性校驗

2.1 研究對象

本文以某萬噸級油船螺旋槳模型為研究對象，該槳模已在上海船舶運輸科學研究所(SSSRI)空泡水筒做過大量試驗[12].螺旋槳模型直徑(D)為248 mm，盤面比為0.4，螺距比(0.35D)為0.826，槳葉數(shù)為4，旋向為右旋.

2.2 網(wǎng)格劃分及可信性校驗

計算區(qū)域與試驗條件完全相同，SSSRI空泡水筒工作段長2.6 m，橫截面呈正方形帶0.1 m倒圓角，寬為0.6 m，高為0.6 m，螺旋槳位于距入口0.5 m處，見圖1，原點坐標位于入口軸系中心處.

圖1 計算區(qū)域(m)Fig.1 Computational field (m)

網(wǎng)格采用六面體核心網(wǎng)格Trim類型，該模型內(nèi)置有一個用以控制表面模板單元格尺寸的上限，這個上限可避免嚴重的模型幾何失真和生成質(zhì)量差的網(wǎng)格.整個計算域分為靜止區(qū)域和螺旋槳旋轉(zhuǎn)區(qū)域，兩個區(qū)域通過交界面連接.為能捕捉到梢渦空泡，需對槳葉梢附近(渦量較大的區(qū)域)進行網(wǎng)格加密.渦量大小的判斷采用Q準則：

(10)

圖2 Q等值面(Q=5×104 s-2)Fig.2 Isosurface of Q(Q=5×104 s-2)

式中：Ω為渦量；S為漩渦變形率，均為二階張量.圖2給出了Q等值面圖，從圖中可看出，螺旋渦管從槳葉導邊起始，由葉梢處脫出，與試驗梢渦空泡起始的位置相同.再進行空化流場計算時發(fā)現(xiàn)，當葉背片空泡產(chǎn)生時，由于片空泡具有一定的厚度，此時渦管在槳葉上的脫出位置與無空化狀態(tài)有所改變，產(chǎn)生向x軸負方向的微小位移.基于上述現(xiàn)象在對該渦管進行網(wǎng)格加密時，加密區(qū)域應包圍此渦管，為降低網(wǎng)格數(shù)量，僅對旋轉(zhuǎn)區(qū)域內(nèi)螺旋渦管進行了加密.

(1)靜止區(qū)域最大網(wǎng)格單元尺寸為0.08D，最小尺寸為0.008D(交界面網(wǎng)格尺寸).為保證槳盤面處伴流場，從入口至下游0.7 m處，以軸系為中心1.62D直徑范圍內(nèi)進行網(wǎng)格加密，網(wǎng)格尺寸為0.02D，網(wǎng)格單元總數(shù)約為123 萬.

圖3 計算域網(wǎng)格(網(wǎng)格2)Fig.3 Mesh of computational field (mesh 2)

(2)旋轉(zhuǎn)區(qū)域最大網(wǎng)格單元尺寸為0.008D，最小網(wǎng)格尺寸為0.002D(槳葉表面)，槳葉設(shè)置為無滑移壁面，邊界層內(nèi)第1層網(wǎng)格高度為2.76×10-4D，共6層，總厚度為4.03×10-4D，本文計算工況下漿葉表面無量綱法向距離y+均在30～100的范圍內(nèi).為捕捉梢渦空化，對槳葉梢及螺旋管進行網(wǎng)格加密，尺寸為0.001D，網(wǎng)格單元總數(shù)約為277 萬.網(wǎng)格劃分見圖3.

圖4 試驗測得的軸向伴流分布Fig.4 Distribution of axial wake measured by the test

為對計算結(jié)果進行可信性校驗，本文運用安全系數(shù)法[13]對3套網(wǎng)格計算得到的螺旋槳平均推力系數(shù)KT進行不確定度估算，KT=T/(ρln2D4)(T為螺旋槳推力，n為螺旋槳轉(zhuǎn)速).網(wǎng)格1、2、3計算得到的KT分別為 0.139 3、0.138 8、0.138 5，試驗值為 0.138 0，不確定度估算結(jié)果為4.464×10-4.分析計算結(jié)果可知，隨著網(wǎng)格數(shù)量的增加，推力系數(shù)逐漸降低且趨于試驗值，但是網(wǎng)格1與2、網(wǎng)格2與3以及網(wǎng)格3與試驗值之間的差異在變小.不確定度值相當小，可認為計算結(jié)果可信.圖5給出了不同網(wǎng)格計算得到的空泡形態(tài).分析圖5可知網(wǎng)格1不能有效地捕捉梢渦空泡，網(wǎng)格2和3均能捕捉到梢渦空泡，葉背片空泡正投影面積和梢渦空泡形態(tài)一致.綜合考慮計算消耗和模擬精度之后，選用網(wǎng)格2進行后續(xù)計算.

圖5 不同網(wǎng)格模擬的空泡形態(tài)與試驗結(jié)果對比Fig.5 Comparisons of cavity pattern between calculation and test results at different mesh sizes

3 結(jié)果分析

3.1 無空化流場數(shù)值模擬

螺旋槳數(shù)值計算條件、工況與試驗條件完全相同[12].試驗中n為30 r/s，KT為0.138，背景壓力為pref，通過改變背景壓力滿足轉(zhuǎn)速空泡數(shù)σn=(pref-p0)/(0.5ρln2D2)=2.96.在正式計算之前，無螺旋槳情況下，需確認在速度進口處給定的流速分布，流向槳盤面位置時速度分布是否發(fā)生改變.圖6給出了入口處和槳盤面處的軸向伴流分布，分析圖6可知在槳盤面處流場分布與入口幾乎相同，因此可按該網(wǎng)格劃分方法進行后續(xù)計算.

為保證推力系數(shù)與試驗相同，又能捕捉到梢渦空化，因此先進行無空化流場計算.采用滑移網(wǎng)格方法，不加入空化模型，計算時間步長設(shè)定為 2.777 8×10-4s，對應螺旋槳旋轉(zhuǎn)角度為3°.背景壓力設(shè)置為1個大氣壓，調(diào)整來流速度(軸向伴流分布保持不變)直至時間平均推力系數(shù)相等.無空化狀態(tài)KT為 0.138 8，誤差為0.6%.螺旋槳旋轉(zhuǎn)一圈推力系數(shù)變化見圖7，分析圖7可知一個旋轉(zhuǎn)周期內(nèi)推力系數(shù)出現(xiàn)4個峰值，與槳葉數(shù)相對應.由圖3可知從θ=45° 附近開始流場逐漸趨于均勻，推力系數(shù)曲線也在此出現(xiàn)波谷，無空化狀態(tài)波谷θ分別為40°、130°、220°及310°.

圖6 不同斷面處軸向伴流分布Fig.6 Distribution of axial wake at different cross-sections

圖7 螺旋槳旋轉(zhuǎn)一圈KT變化Fig.7 Variations of KT during one propeller rotation

3.2 空化流場數(shù)值模擬

當無空化流場計算穩(wěn)定后，激活歐拉多相流模型和空化模型，調(diào)整環(huán)境壓力使轉(zhuǎn)速空泡數(shù)為2.96，計算時間步長為 9.259 26×10-5s，對應螺旋槳旋轉(zhuǎn)角度為1°.空化狀態(tài)KT為0.139，誤差為0.72%，推力脈動高于無空化狀態(tài)，最小推力系數(shù)出現(xiàn)在θ=45°，135°，225°，315°.為驗證數(shù)值計算結(jié)果的準確性，與試驗結(jié)果比較，對不同相位角槳葉空泡形態(tài)進行了對比分析，見圖8.采用αv=0.1等值面描述空泡形態(tài).圖中給出了槳葉進、出伴流區(qū)域空泡發(fā)展及潰滅的整個過程：當槳葉開始進入伴流區(qū)域時首先出現(xiàn)葉背片空泡，起始于槳葉導邊附近，隨著槳葉繼續(xù)旋轉(zhuǎn)，葉背片空泡面積逐漸增大，梢渦空泡開始脫出，最大葉背片空泡面積出現(xiàn)在θ=10° 的位置，隨后片空泡面積逐漸減小，梢渦空化逐漸加強，隨著槳葉逐漸出伴流區(qū)域，葉背片空泡消失，梢渦逐漸降弱直至消失.數(shù)值計算準確的預測了這一歷程，片空泡初生、片空泡面積、梢渦初生及梢渦消失等現(xiàn)象與試驗結(jié)果吻合.

但是對于梢渦空化而言，空化形態(tài)略有差異，例如θ=60° 時，試驗照片顯示槳葉脫出的梢渦并不平滑，剛從葉梢梢部脫落時，梢渦非常細，隨后梢渦空化加強，渦管變粗且外圍呈現(xiàn)霧狀空化，渦管表面非常粗糙.這可能是由于梢渦脫出后進入高伴流區(qū)域，速度梯度變大，壓力非常不穩(wěn)定，致使空化呈現(xiàn)出強烈的非定常特性，在短時間內(nèi)劇烈變化.計算中θ=60° 時梢渦空化沿螺旋渦管也是由弱變強，從葉梢附近脫出、消失再出現(xiàn)，但是梢渦空間結(jié)構(gòu)更加平順光滑.這主要是由于網(wǎng)格尺寸不能無限增加，再者RANS方法對該微小尺度渦不能夠準確預測.即使采用LES方法，其梢渦渦管依然平順[7].

試驗中空泡手描圖能夠?qū)θ~背片空泡的大小、形態(tài)及穩(wěn)定性進行詳盡描述，從θ=20° 開始梢渦在空間內(nèi)呈現(xiàn)復雜的三維結(jié)構(gòu)，手描圖難以在二維平面內(nèi)對其準確的描述，對于梢渦形態(tài)手描圖僅進行了示意，因此僅對葉背片空泡進行了量化計算，分析了與試驗結(jié)果的差異性.現(xiàn)將不同相位角空泡計算結(jié)果、試驗手描圖和試驗照片進行圖像處理，每張照片均進行圖像重塑為正投影圖片，分別測繪葉背空泡正投影平面面積和槳葉正投影平面面積.定義葉背片空泡面積百分比γ=AS/AB(AS為葉背空泡面積，AB為槳葉面積).表1給出了數(shù)值和試驗葉背片

表1 計算和試驗葉背片空泡面積對比Tab.1 Comparisons of back-sheet cavity area between calculation and test results

圖8 不同θ時空泡形態(tài)計算結(jié)果與試驗結(jié)果比較Fig.8 Comparisons of cavity pattern between calculation results and test results observation at different θ values

空泡正投影面積百分比結(jié)果，分析表1可知數(shù)值計算結(jié)果與試驗結(jié)果吻合，僅在θ=40° 時相差5%左右，其他情況下計算結(jié)果與試驗結(jié)果的差值在2.5%以內(nèi).

4 結(jié)語

本文以較少的網(wǎng)格數(shù)量成功模擬了非均勻流場螺旋槳葉背片空泡和梢渦空泡.結(jié)果表明：Schnerr-Sauer空化模型和可實現(xiàn)的k-ε兩層湍流模型能夠?qū)Ψ蔷鶆蛄鲌雎菪龢张菪螒B(tài)進行準確預測，計算結(jié)果與試驗結(jié)果吻合，該計算方法可為非均勻流場螺旋槳空泡模擬提供參考.