非均勻流中E779A螺旋槳空化渦流演變及壓力脈動數(shù)值分析

何朝暉，龍 云，韓承灶，季 斌

（武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072）

0 引言

近年來，我國艦船往大型化和高速化發(fā)展，這將導(dǎo)致船用螺旋槳負(fù)荷的不斷增加，從而產(chǎn)生空化?？栈粌H會降低螺旋槳的水力性能，而且會加劇流場的不穩(wěn)定性，誘導(dǎo)周期性的壓力脈動，引起噪聲以及空蝕。通過數(shù)值模擬方法研究空化對流場的影響對于改善船舶性能有著指導(dǎo)作用。關(guān)于空化流動的數(shù)值模擬，研究人員最開始研究的是幾何形狀簡單的翼型和回轉(zhuǎn)體，早期螺旋槳的空化流動模擬研究采用勢流方法，如升力面方法、面元法等進(jìn)行求解，隨著粘流CFD技術(shù)的發(fā)展，到20世紀(jì)90年代，螺旋槳的空化流動CFD 數(shù)值模擬研究逐漸出現(xiàn)并得到快速發(fā)展。與勢流理論求解不同，粘流CFD 數(shù)值模擬方法將粘性的影響加以考慮，能夠更為深入地對空化進(jìn)行描述。Bensow等[1]通過自定義公式近似非均勻來流速度分布，模擬研究了E779A 螺旋槳的非定?？栈鲃樱籗hin等[2]基于HEM 方法和蒸汽輸運(yùn)方程，模擬了船后常規(guī)及大側(cè)斜螺旋槳周圍的空化流動，結(jié)果表明非定?？栈螒B(tài)在定性上精度可接受；朱志峰等[3]采用了混合多相流模型，求解了均勻入流工況下的E779A螺旋槳的空化性能，較為準(zhǔn)確地預(yù)報了片空化；葉金銘等[4]基于空泡混合模型，采用RANS 方法模擬計算了非均勻來流的螺旋槳的空化性能，空泡計算結(jié)果與試驗較為吻合。近年來盡管數(shù)值模擬方法有所改進(jìn)，但許多研究仍然局限于螺旋槳性能和空化形態(tài)的預(yù)測。本文以E779A 螺旋槳為研究對象，該槳于意大利船模水池中心經(jīng)過了詳細(xì)的敞水性能和空化特性的試驗研究。此外，在歐盟相關(guān)項目中曾組織多家單位對其均勻來流空化進(jìn)行了數(shù)值模擬比較，擁有較為完整詳細(xì)的數(shù)據(jù)。本文重點關(guān)注了螺旋槳空化流場的壓力脈動及渦量演化，研究了空化對壓力脈動以及空化對渦流場的影響。

1 數(shù)值模擬方法

本文采用基于求解RANS 方程的方法，結(jié)合k-ωSST 湍流模型和Zwart-Gerber-Belamri 空化模型開展螺旋槳空化流動數(shù)值模擬研究。

1.1 湍流模型

k-ωSST 湍流模型由Menter 提出，用于改進(jìn)對逆壓梯度流的預(yù)測，其能夠較好地預(yù)測旋轉(zhuǎn)機(jī)械周圍的空化流動[5]。k-ωSST湍流模型使用混合函數(shù)F1將近壁區(qū)的k-ω模型與遠(yuǎn)場的k-ε模型混合。湍流動能k和湍流頻率ω的傳輸方程為

式中參數(shù)計算和常數(shù)的選取可參考相關(guān)文獻(xiàn)。

1.2 空化模型

本文采用Zwart-Gerber-Belamri（以下簡稱ZGB）空化模型[6]來捕捉由空化引起的兩相流特征，質(zhì)量傳遞用右側(cè)帶有兩個源項的蒸汽傳輸方程來處理

單個氣泡中的質(zhì)量變化率與簡化的Rayleigh-Plesset方程相結(jié)合，以閉合方程(3)，蒸發(fā)和冷凝源項最終表示為

式中：ρv代表氣體密度，α代表氣體體積分?jǐn)?shù)，t代表時間，u代表速度，αnuc代表成核區(qū)的體積分?jǐn)?shù)，Rb代表氣泡半徑，ρl代表液體密度，F(xiàn)vap和Fcond分別為ZGB空化模型的經(jīng)驗蒸發(fā)系數(shù)和凝結(jié)系數(shù)。常用的參數(shù)取值有：Rb=1μm，αnuc=5E-4，F(xiàn)vap=50，F(xiàn)cond=0.01。

2 幾何、網(wǎng)格和邊界條件

E779A 螺旋槳為四葉槳，槳模直徑為0.227 27 m，其計算域和邊界條件見圖1。螺旋槳位于旋轉(zhuǎn)域，槳葉表面設(shè)置為無滑移壁面。靜止部分的外圓柱表面設(shè)置為自由滑移壁面。入口邊界的非均勻來流速度借助CFX 求解器里的“profile”功能導(dǎo)入試驗數(shù)據(jù)獲得，見圖2。螺旋槳的轉(zhuǎn)速設(shè)置為30.5 r/s。螺旋槳的進(jìn)速系數(shù)為

圖1 非均勻來流計算域和邊界條件Fig.1 Computational domain and boundary condition of non-uniform inlet flow

圖2 數(shù)值模擬中的非均勻來流速度分布Fig.2 Non-uniform inflow velocity distribution in the simulation

式中，n為螺旋槳轉(zhuǎn)速，D為螺旋槳直徑。

非均勻進(jìn)口邊界與螺旋槳之間的距離設(shè)置為0.7D，出口為螺旋槳下游6D處，設(shè)置為靜壓出口，其壓力由空化數(shù)求得，公式如下：

式中，Pv表示水的飽和蒸汽壓力，ρ代表水的密度，pout表示流域出口壓力。

通過ANSYS ICEM CFD 軟件采用相同的拓?fù)浣Y(jié)構(gòu)生成了網(wǎng)格數(shù)逐漸增大的三套結(jié)構(gòu)化網(wǎng)格，三套網(wǎng)格在x、y、z方向均以 2 倍進(jìn)行遞增，網(wǎng)格總數(shù)分別為186萬、502萬和1404萬。選取非均勻來流、進(jìn)速系數(shù)為0.9、空化數(shù)為4.455 的空化工況，對E779A 螺旋槳的空化流動進(jìn)行了數(shù)值模擬，使用三套網(wǎng)格計算得到的結(jié)果見圖3。對于圖中選定的三個典型時刻，當(dāng)網(wǎng)格數(shù)量較少時（網(wǎng)格1），數(shù)值模擬得到的空化形態(tài)與文獻(xiàn)[7]中試驗結(jié)果相比，雖然在變化趨勢上是吻合的，但是空化區(qū)域的長度和大小均偏小。而隨著網(wǎng)格的逐漸細(xì)化，預(yù)測的空化面積明顯增加，與試驗結(jié)果更加吻合，但當(dāng)網(wǎng)格細(xì)化到一定程度后，改善效果不明顯。因此，在保持計算精度的前提下，為了節(jié)約計算資源，本文所有的計算選取網(wǎng)格2進(jìn)行計算。

圖3 對比3套網(wǎng)格預(yù)測的空化形態(tài)Fig.3 Comparison of cavitation patterns predicted by three sets of grids

E779A 螺旋槳及其計算域的結(jié)構(gòu)化網(wǎng)格見圖4～5。為了捕獲相關(guān)位置詳細(xì)的流場信息，在導(dǎo)邊、隨邊及葉梢區(qū)域進(jìn)行了局部加密。無空化流動的收斂結(jié)果被設(shè)定為非定?？栈鲃佑嬎愕某跏剂鲌?。無空化流動殘差收斂標(biāo)準(zhǔn)設(shè)置為10-5。非定常計算內(nèi)部迭代的最大次數(shù)為每個時間步長20 次，收斂標(biāo)準(zhǔn)采用10-4來平衡時間消耗和計算精度。

圖4 E779A螺旋槳結(jié)構(gòu)化網(wǎng)格Fig.4 Structured meshes of the E779A propeller

圖5 E779A螺旋槳非均勻來流計算域結(jié)構(gòu)化網(wǎng)格Fig.5 Structured meshes of the fluid domain of non-uniform inlet flow

3 結(jié)果與討論

3.1 非定?？栈螒B(tài)

由于船體幾何的干擾，螺旋槳前方水流的來流速度分布往往是不均勻的。為了使數(shù)值模擬的結(jié)果貼合實際，本文對E779A 螺旋槳在非均勻來流下的空化流動進(jìn)行了數(shù)值模擬，參照文獻(xiàn)[7]的計算工況，設(shè)置進(jìn)速系數(shù)為0.9，空化數(shù)為4.455。由于葉片載荷的周期性變化，在E779A 螺旋槳葉片上觀察到空化形態(tài)的瞬態(tài)變化，見圖6，圖中展示了E779A 螺旋槳在數(shù)值模擬中不同時刻的汽相體積分?jǐn)?shù)為0.1 的等值面圖，并與文獻(xiàn)[7]中給出的試驗得到的空化演變圖進(jìn)行比較。

圖6 右圖中展示了E779A 螺旋槳在試驗中的空化演變圖。從圖中可以看到：對于選定的槳葉，在旋轉(zhuǎn)角為-25°時，螺旋槳葉背的導(dǎo)邊近葉梢處，空化得到初步的形成，并隨著角度的增大，空化的體積在逐漸增大；在旋轉(zhuǎn)角為0°時，空化的發(fā)展幾乎達(dá)到最大程度，此時，片空化逐漸卷起進(jìn)入梢渦，空化開始變得不穩(wěn)定；在旋轉(zhuǎn)角為10°時，葉片已經(jīng)開始離開尾流速度較小的區(qū)域，空化從下邊區(qū)域開始收縮潰滅，體積有所減小，并且可以觀察到延伸至梢渦區(qū)域的空化旋轉(zhuǎn)運(yùn)動；在旋轉(zhuǎn)角為20°時，空化進(jìn)一步往葉梢區(qū)域收縮，大部分都轉(zhuǎn)變成了不穩(wěn)定的云空化，梢渦空化得到進(jìn)一步發(fā)展。

圖6 E779A螺旋槳非均勻來流工況數(shù)值模擬與文獻(xiàn)[7]空化形態(tài)對比Fig.6 Comparison of the predicted and experimental cavity patterns of E779A propeller of non-uniform inlet flow

從圖6 左圖中可以看到：對于選定的槳葉，在旋轉(zhuǎn)角為-25°時，空化初步在葉背導(dǎo)邊處形成，并隨著角度的增大，空化的體積逐步增大；在旋轉(zhuǎn)角為0°時，空化的體積達(dá)到最大值，此后，隨著角度的增大，空化的體積逐漸向上縮?。辉谛D(zhuǎn)角度為20°時，可以明顯觀察到梢渦空化的形成。可見，本文對E779A 螺旋槳的非定常模擬結(jié)果與試驗結(jié)果相比，無論在空化的形態(tài)還是空化的演變趨勢上均有較好的一致性，這表明本文的數(shù)值模擬結(jié)果是可靠的。

3.2 壓力脈動

上文對非定?？栈螒B(tài)的分析，驗證了本文非定常空化數(shù)值模擬方法的可靠性。為了進(jìn)一步研究空化的產(chǎn)生對流場的影響，依據(jù)試驗[8]中的布置，本文選取距離槳轂中心300 mm 的四個監(jiān)測點對螺旋槳的壓力脈動進(jìn)行了分析，其空間位置見圖7。數(shù)值模擬得到的螺旋槳模型在無空化與空化條件下，四個不同位置監(jiān)測點（p1、p2、p3、p4）處的壓力脈動，見圖8。圖中，壓力脈動幅值表示為無量綱參數(shù)Kp[9]：

圖7 E779A螺旋槳壓力脈動監(jiān)測點示意圖Fig.7 Sketch of pressure fluctuation monitoring points for E779A propeller

式中，p為監(jiān)測點壓力，ρ為流體密度，Kb為固壁參數(shù)，取為2。

圖8（a）為無空化工況，四個監(jiān)測點處的脈動壓力頻率是相同的，均為螺旋槳的葉片通過頻率。脈動壓力系數(shù)在4.552～4.582的范圍內(nèi)上下波動，波動范圍較小。對比四幅曲線圖可以看到，四個監(jiān)測點具有的脈動壓力相位是不相同的，這是因為各監(jiān)測點相對于螺旋槳的幾何位置不同。圖8（b）為空化工況，脈動壓力系數(shù)波動范圍增大，在4.51～4.61之間。對比四幅圖可以看到，空化發(fā)生后，不僅四個監(jiān)測點脈動壓力的變化頻率同為葉片的通過頻率，而且其壓力脈動的相位也是一樣的。這正是空化條件下的脈動壓力變化趨勢與無空化條件下的脈動壓力變化趨勢最明顯的區(qū)別。

圖8 四個監(jiān)測點上的壓力脈動Fig.8 Pressure fluctuations at four monitoring points

計算得到的四個監(jiān)測點的壓力脈動幅值變化曲線，見圖9。在無空化工況（圖9（a）），四個監(jiān)測點處的壓力脈動主要成分是一階葉片通過頻率分量，脈動幅值均較小，在0.004 以下。而在空化發(fā)生后（圖9（b）），壓力脈動主要成分除了一階葉片通過頻率分量，還有二階乃至更高階葉片通過頻率分量。通過圖9（a）和圖9（b）的對比可以明顯看到，空化發(fā)生后，四個監(jiān)測點的一階葉片通過頻率分量顯著增大，均在0.01以上。而且由于空化的影響，二階葉片通過頻率分量也遠(yuǎn)比無空化時的一階葉片通過頻率大得多，這表明空化時的壓力脈動幅值遠(yuǎn)強(qiáng)于無空化時的壓力脈動幅值。此外，由于螺旋槳的每個槳葉都會經(jīng)歷空化的周期性生長和潰滅的演變過程，這個過程會使流場中誘發(fā)大幅度的壓力脈動，導(dǎo)致所有監(jiān)測點處壓力脈動的同步增大。因此，一旦船后的螺旋槳上發(fā)生了周期性演變的空化，就會引發(fā)對船體有害的振動和噪聲，從而影響船舶運(yùn)行的舒適性和安全性。

圖9 四個監(jiān)測點上的壓力脈動幅值Fig.9 Amplitude of pressure fluctuation at four monitoring points

為了定量評估非定?？栈鲃酉聰?shù)值計算得到的壓力脈動的結(jié)果，將四個監(jiān)測點的壓力脈動幅值的一階分量與相對應(yīng)的試驗結(jié)果[8]進(jìn)行分析比較，見圖10。圖中數(shù)值結(jié)果及試驗結(jié)果通過公式（8）進(jìn)行了無量綱化。從圖中可以看到，數(shù)值模擬的結(jié)果均比試驗的結(jié)果要稍稍偏大，兩者的平均誤差為8.93%，在p4 點誤差最大，為16.89%，但總體來說計算結(jié)果與試驗結(jié)果是較為吻合的。這些結(jié)果表明，使用當(dāng)前的數(shù)值方法預(yù)測螺旋槳周圍場點的壓力脈動是可行的。這有助于工程上預(yù)防及控制空化引起的壓力脈動帶來的不利影響。

圖10 四個監(jiān)測點壓力脈動幅值的一階葉片通過頻率分量的對比Fig.10 Comparison of the first blade passing frequency component of pressure amplitude at four monitoring points

3.3 渦量輸運(yùn)方程

為了進(jìn)一步研究空化對螺旋槳周圍渦流場的影響，采用圓柱坐標(biāo)系[10]中的相對渦量輸運(yùn)方程進(jìn)行分析，其表達(dá)式如下：

式中，ω是相對渦量，u是相對速度，v是運(yùn)動粘度，ρ是混合密度，ω1是螺旋槳的旋轉(zhuǎn)速度。

在此方程中，渦量變化率由右側(cè)的五個分量所組成。其中第一項代表相對渦流拉伸，由速度梯度決定；第二項是由于體積變化引起的相對渦流收縮膨脹，這說明了流體可壓縮性對渦流的影響；第三項是斜壓扭矩項，由壓力梯度與密度梯度不平衡而產(chǎn)生，在正壓流體中等于零；第四項是粘性擴(kuò)散項，反映流體相對運(yùn)動中因粘性而導(dǎo)致的渦量擴(kuò)散；第五項表明了科氏力在旋轉(zhuǎn)坐標(biāo)系中對相對渦量的影響。

空化氣相體積分?jǐn)?shù)αv=0.1 等值面圖、氣相體積分?jǐn)?shù)αv=0.1 的云圖及上式中渦量輸運(yùn)方程右側(cè)五項周向方向分量的云圖見圖11。圖中沿空化的發(fā)展方向選取了6個觀察平面，分別位于空化的前部，中部以及尾部。從圖中可以看到，渦量場的演變與空化是密切相關(guān)的。相對周向渦量拉伸項與其他四項相比明顯是較大的，而隨著空化的發(fā)展，葉尖區(qū)域的相對渦量拉伸項呈現(xiàn)出更加復(fù)雜的發(fā)展趨勢，這表明空化的發(fā)展會促進(jìn)拉伸項的發(fā)展，從而加劇湍流和渦流的不穩(wěn)定性。周向相對渦量拉伸項在空化內(nèi)外的分布呈現(xiàn)正負(fù)交錯，這表明沿該區(qū)域相對較強(qiáng)的渦量輸運(yùn)對旋渦的發(fā)展具有促進(jìn)作用。在螺旋槳的尾流附近，周向相對渦量拉伸項仍然很大，是影響渦量輸送的主要因素。此外，與速度散度相關(guān)的收縮膨脹項主要集中在空化內(nèi)部，表明了空化內(nèi)部流動可壓縮性的影響。圖中清楚地顯示了周向相對渦斜壓扭矩項沿汽液交界面是顯著的，這是由該處密度梯度和壓力梯度不平衡所導(dǎo)致的。在非定常流動中，由空化引起的速度散度和壓力、密度梯度不平衡使得收縮膨脹項和斜壓扭矩項都是促進(jìn)渦量產(chǎn)生和改變渦量場的重要源項。而剩余的相對軸向渦粘性耗散項和科氏力項相對于其他三項而言比較小，對渦流場發(fā)展的貢獻(xiàn)較小。綜上所述，空化發(fā)生后，引起的渦量拉伸、收縮膨脹及斜壓力矩會促進(jìn)渦的產(chǎn)生和發(fā)展，加劇流場的不穩(wěn)定。

圖11 預(yù)測的相對渦量拉伸項，膨脹項，斜壓扭矩項，粘性擴(kuò)散項，科氏力項Fig.11 Predicted relative vortex stretching,vortex dilatation,baroclinic torque and Coriolis force terms

4 結(jié) 論

本文采用求解基于k-ωSST 湍流模型和Zwart-Gerber-Belamri 空化模型的RANS 方程的方法對E779A螺旋槳的非定?？栈鲃舆M(jìn)行了數(shù)值模擬，主要結(jié)論如下：

（1）通過對E779A 螺旋槳的非均勻來流空化進(jìn)行模擬研究，得到的非定?？栈療o論形態(tài)還是演變過程均與試驗結(jié)果較為吻合，驗證了本文計算的可靠性。

（2）對螺旋槳周圍無空化及空化壓力脈動的對比分析結(jié)果表明，空化的發(fā)生將會誘發(fā)大幅度周期性演變的壓力脈動，從而引發(fā)對船舶不利的振動和噪聲。

（3）對空化流場中相對渦量輸運(yùn)方程各項分量的比較分析結(jié)果表明，由空化引起的渦量拉伸、收縮膨脹和斜壓力矩是影響渦量場的主要因素，其會促進(jìn)渦的生成和發(fā)展，從而加劇流場的不穩(wěn)定性。