劉 恒，伍 銳,3，孫 碩

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

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

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

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

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

1 控制方程和空泡模型

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

(1)

(2)

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

(3)

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

(4)

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

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

(5)

(6)

(7)

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

當(dāng)p

(8)

當(dāng)p>p0時(shí)，氣體凝結(jié)，質(zhì)量輸運(yùn)速率為

(9)

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

2.1 研究對(duì)象

本文以某萬(wàn)噸級(jí)油船螺旋槳模型為研究對(duì)象，該槳模已在上海船舶運(yùn)輸科學(xué)研究所(SSSRI)空泡水筒做過(guò)大量試驗(yàn)[12].螺旋槳模型直徑(D)為248 mm，盤(pán)面比為0.4，螺距比(0.35D)為0.826，槳葉數(shù)為4，旋向?yàn)橛倚?

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

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

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

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

(10)

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

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

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

圖3 計(jì)算域網(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è)置為無(wú)滑移壁面，邊界層內(nèi)第1層網(wǎng)格高度為2.76×10-4D，共6層，總厚度為4.03×10-4D，本文計(jì)算工況下漿葉表面無(wú)量綱法向距離y+均在30～100的范圍內(nèi).為捕捉梢渦空化，對(duì)槳葉梢及螺旋管進(jìn)行網(wǎng)格加密，尺寸為0.001D，網(wǎng)格單元總數(shù)約為277 萬(wàn).網(wǎng)格劃分見(jiàn)圖3.

圖4 試驗(yàn)測(cè)得的軸向伴流分布Fig.4 Distribution of axial wake measured by the test

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

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

3 結(jié)果分析

3.1 無(wú)空化流場(chǎng)數(shù)值模擬

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

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

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

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

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

表1 計(jì)算和試驗(yàn)葉背片空泡面積對(duì)比Tab.1 Comparisons of back-sheet cavity area between calculation and test results

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

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

4 結(jié)語(yǔ)

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