基于面元-渦粒子法的螺旋槳氣動特性及噪聲研究

劉乾, 劉漢儒, 李家輝, 尚珣, 陳南樹, 王掩剛

1.西北工業(yè)大學(xué) 動力與能源學(xué)院, 陜西 西安 710072;2.中國空氣動力研究發(fā)展中心 氣動噪聲控制重點實驗室, 四川 綿陽 621000

隨著航空業(yè)的發(fā)展,對飛行器的高效、節(jié)能、低噪聲要求越來越高。渦扇發(fā)動機(jī)對燃料的利用效率約40%,而分布式電推進(jìn)系統(tǒng)對電能的利用率能夠超過70%[1-2],小型螺旋槳在低速設(shè)計工況下能長期保持較高氣動效率。要使分布式推進(jìn)系統(tǒng)的噪聲滿足標(biāo)準(zhǔn),需要使用非定常氣動及聲學(xué)仿真手段進(jìn)行預(yù)測,并采取降噪措施處理使之滿足適航規(guī)定。螺旋槳氣動噪聲主要來源于空間中的非定常壓力脈動和雷諾應(yīng)力[3],為了實現(xiàn)葉輪機(jī)械噪聲性能的預(yù)測、評估并提出改進(jìn)的方向,必須對葉輪機(jī)械開展非定常數(shù)值計算。

面元法是目前計算精度最高的勢流方法[4-6],Baltazar等[7]使用面元法對導(dǎo)管螺旋槳性能進(jìn)行計算,發(fā)現(xiàn)尾渦計算準(zhǔn)確性對性能預(yù)測極為重要。為了準(zhǔn)確模擬尾跡作用,渦粒子法被研究人員用作尾跡模型。渦粒子法采用離散的拉格朗日粒子求解渦量輸運(yùn)方程[8-9],數(shù)值耗散更低,能夠準(zhǔn)確模擬尾跡的發(fā)展情況[10],將面元法和渦粒子法結(jié)合可以實現(xiàn)計算速度與精度的統(tǒng)一。Willis等[11]使用面元-渦粒子法對翼身振動問題進(jìn)行了求解,并且獲得了準(zhǔn)確的結(jié)果。胡昊等[12]使用面元-渦粒子法對風(fēng)力機(jī)氣動特性進(jìn)行計算,發(fā)現(xiàn)面元-渦粒子法能夠給出風(fēng)力機(jī)的其他流動細(xì)節(jié)。

在上述研究中,將面元法與渦粒子法結(jié)合進(jìn)行葉輪機(jī)械氣動性能模擬被認(rèn)為是具有研究價值的工作,但將該方法用于氣動噪聲預(yù)測的研究有限。面元法高效的非定常數(shù)值模擬性能,使之在聲學(xué)預(yù)測中具備獨特優(yōu)勢,把面元-渦粒子法與氣動聲學(xué)模型結(jié)合,對葉輪機(jī)械氣動及聲學(xué)性能的預(yù)測有較高的工程實用價值,適于葉輪機(jī)械的多工況、多學(xué)科設(shè)計。

本文將分別介紹面元法、渦粒子法、Lowson頻域聲學(xué)模型以及面元-渦粒子法的理論,結(jié)合有限體積法驗證面元-渦粒子法的仿真精度,并深入對比分析2種方法的數(shù)值計算結(jié)果,進(jìn)一步揭示螺旋槳氣動特性和聲學(xué)特性。

1 理論方法

1.1 面元法

面元法以流體連續(xù)性為理論基礎(chǔ),求解三維流場中勢函數(shù):

2φ=0

(1)

使用Green函數(shù)求解(1)式,并引入源、匯、偶極等作為勢流單元和不可穿透固體表面S,則在S上的P,Q兩點間存在:

(2)

根據(jù)固體壁面、尾跡面元不可穿透邊界條件,無限遠(yuǎn)處勢函數(shù)誘導(dǎo)速度為0假設(shè)

(3)

可使用偶極和源表示速度勢,從(3)式得到:

本文中使用渦線段等效面元計算誘導(dǎo)速度。如果渦線段與場點較近會產(chǎn)生數(shù)值計算的奇異性,因此引入渦環(huán)半徑r減弱奇異性

式中:R1是渦線段起點和場點間的空間向量;R2是渦線段終點和場點間的空間向量;R0是渦線段起點和終點間的空間向量;h是渦線段和場點間的垂直距離。從(6)式能看出,當(dāng)h過小時,r的存在將會減弱誘導(dǎo)速度數(shù)值奇異性;計算較遠(yuǎn)處誘導(dǎo)速度時,r/h接近0,對計算精度基本無影響。

為了保證面元能夠準(zhǔn)確捕捉前緣、尾緣處的流動情況,應(yīng)該有效加密面元。為了更好地滿足Kutta條件,把圓尾緣處理為尖尾緣,并使用3次樣條插值方法生成離散坐標(biāo)點。

本文使用余弦方法對螺旋槳的徑向進(jìn)行面元劃分,控制吸、壓力面的周向和徑向離散節(jié)點數(shù)量一致,使用余弦方法離散可以保證前、尾緣處的面元密度更大,能準(zhǔn)確描述葉片曲率。螺旋槳吸力面或壓力面的徑向、弦向面元數(shù)分別為NR和NC,所以單個螺旋槳葉片可被劃分為2×NR×NC個面元。以NACA0012二維標(biāo)準(zhǔn)翼型示意,使用余弦方法劃分面元[13]。

圖1 NACA0012的余弦離散

1.2 渦粒子法

渦粒子法使用離散的拉格朗日粒子求解渦量輸運(yùn)方程。對不可壓Navier-Stokes方程取散度后可以得到渦量輸運(yùn)方程

(7)

使用渦粒子離散空間中的渦量,設(shè)第i個渦粒子的強(qiáng)度和位置分別為αi和xi,則渦量場離散為

(8)

式中:渦粒子強(qiáng)度αi是渦量和體積的乘積;ξσ是使渦量光滑分布的磨光算子(光滑核函數(shù)),普遍使用高斯形式的核函數(shù)(也稱截斷函數(shù))[14]

(9)

使用渦粒子模擬流場的發(fā)展時,粒子強(qiáng)度因為旋渦拉伸和黏性耗散發(fā)生變化,因此渦粒子法的求解過程分為位置和渦強(qiáng)兩部分之間的迭代

(10)

(11)

(11)式右邊由旋渦拉伸項和黏性擴(kuò)散項構(gòu)成,分別使用核函數(shù)和渦強(qiáng)交換方法(PSE)對2項求解,得到

使用四階Runge-Kutta算法[15]求解(12)～(13)式。(12)式中ρ=|x-xj|/σ是無量綱長度參數(shù),K(ρ)是用于計算誘導(dǎo)速度的核函數(shù)。

光滑核函數(shù)直接決定空間中渦量的分布情況,所以核函數(shù)的選擇會顯著影響數(shù)值計算精度。根據(jù)矩特性判斷光滑核函數(shù)的精度,對r階精度的光滑核函數(shù)有如下約束[16]:

(16)

給出幾組不同精度的光滑函數(shù)和與之相對應(yīng)的Green函數(shù)，如表1和圖2所示。

表1 常用的光滑函數(shù)和對應(yīng)的Green函數(shù)

圖2 3種不同精度的Green函數(shù)

在表1常用的光滑函數(shù)和對應(yīng)的Green函數(shù)中有

(17)

如圖2所示,隨著階數(shù)的增加,渦量分布越集中,數(shù)值計算的誤差也越小;Fun2和Fun3能把渦量集中在相同范圍內(nèi),但Fun2比Fun3階數(shù)更低,計算量更小,所以使用二階精度的Gaussian分布函數(shù)計算誘導(dǎo)速度。

1.3 Lowson聲學(xué)模型

葉輪機(jī)械中壓力脈動有較強(qiáng)的周期性,使用時域方法計算聲壓要求輸入多個時間步的壓力脈動,計算成本較高;而頻域方法僅需要計算葉片特征頻率諧頻上的聲壓,所以使用頻域下的Lowson聲壓表達(dá)式[17]

cn=an+ibn=

(18)

將時間微分項替換,引入周向和軸向力分量,分部積分后得到Lowson公式

(19)

Fi=(-T,-Dsinθ,Dcosθ)

(xi-yi)=(x,y-Rcosθ,-Rsinθ)

為了提高適用性,將片條理論假設(shè)下的聲壓公式變換到三維坐標(biāo)系[18]

對第n階旋轉(zhuǎn)偶極子聲源的自由遠(yuǎn)場聲壓表達(dá)式有

(22)

1.4 面元-渦粒子法

使用簡化的Hess等效原則,實現(xiàn)面元到渦粒子的轉(zhuǎn)化。為了滿足尾緣的Kutta條件,添加一列緩沖尾跡面元作為面元和渦量的過渡,并將其轉(zhuǎn)化的渦粒子靠近尾緣邊等效為渦線(圖3尾跡面元轉(zhuǎn)化為渦粒子中紅色線段),其余3條邊轉(zhuǎn)化為渦粒子:

圖3 尾跡面元轉(zhuǎn)化為渦粒子

在求解過程中,渦粒子產(chǎn)生的誘導(dǎo)速度作為自由來流速度的一部分,用于求解面元上的源匯分布;面元對各個渦粒子產(chǎn)生誘導(dǎo)速度和速度梯度,作為自由來流一部分參與運(yùn)算。具體迭代過程見圖4。

CFD計算結(jié)果與聲學(xué)模型耦合需要使用非定常數(shù)值計算方法獲得時域下葉表壓力脈動,對時域壓力脈動做時頻變換,將頻域聲壓脈動作為聲源項輸入聲類比模型,流程如圖5所示。

圖4 面元法和渦粒子法的耦合過程

圖5 CFD和聲學(xué)模型耦合

2 數(shù)值計算結(jié)果

2.1 前處理參數(shù)設(shè)置

在軸向來流速度15 m/s、轉(zhuǎn)速為5 000 r/min的工況下(進(jìn)距比為4.643),對螺旋槳進(jìn)行仿真。

2.1.1 有限體積法網(wǎng)格尺度

流場分為旋轉(zhuǎn)域和外域,旋轉(zhuǎn)域半徑取螺旋槳半徑的1.5倍,外域取螺旋槳半徑的10倍。使用URANS和聲類比理論做葉輪機(jī)械氣動噪聲分析時,選擇3BPF作為分析頻率的上限,采樣頻率要高于6BPF使之滿足奈奎斯特采樣定理。

3BPF對應(yīng)的聲波波長是0.243 m,為了能夠嚴(yán)格滿足點聲源假設(shè),網(wǎng)格尺度需要小于1/6聲波波長,所以網(wǎng)格尺度應(yīng)該小于0.040 5 m。以表2中的網(wǎng)格尺度參數(shù)作為基準(zhǔn)。

表2 網(wǎng)格大小設(shè)置參數(shù) m

劃分多面體網(wǎng)格如圖6所示,為了驗證該算例的準(zhǔn)確性,需要進(jìn)行網(wǎng)格無關(guān)性驗證:

1) 邊界層網(wǎng)格高度

在表2所示算例中,控制邊界層厚度為7×10-4m,網(wǎng)格層數(shù)為15層,保證近壁面處網(wǎng)格高度小于10-5m。改變近壁面網(wǎng)格高度后,求得推力相對誤差小于0.1%,說明邊界層網(wǎng)格高度合適。

2) 尾跡區(qū)網(wǎng)格密度

尾跡區(qū)域中,在原網(wǎng)格基礎(chǔ)上疊加5%基礎(chǔ)尺寸網(wǎng)格,求得推力為11.474 N,與基準(zhǔn)算例11.437 N相對誤差為0.32%,認(rèn)為尾跡區(qū)網(wǎng)格密度滿足要求。

圖6 旋轉(zhuǎn)域網(wǎng)格示意圖

分別對邊界層、葉表、尾跡區(qū)進(jìn)行網(wǎng)格加密并求解,證實使用該算例網(wǎng)格可以獲得準(zhǔn)確結(jié)果。

2.1.2 面元離散尺度

結(jié)合圖1所示的面元余弦離散方法,對螺旋槳葉片的吸、壓力面分別進(jìn)行弦向和展向50×20,50×30,50×40,50×50共4種方式離散,并探究不同面元離散密度對數(shù)值仿真的影響。

表3 不同面元離散密度和有限體積法求得推力

圖7 螺旋槳葉片面元離散示意圖

如圖8所示,面元離散密度與計算精度有如下關(guān)系:

1) 隨著展向面元離散密度增加,面元-渦粒子法求出的推力逐漸接近有限體積法計算結(jié)果;

2) 當(dāng)面元密度過高時,面元-渦粒子法的推力偏離有限體積法計算結(jié)果;

3) 面元長寬比接近1時,計算結(jié)果更準(zhǔn)確。

上述現(xiàn)象說明,為了提高計算精度,需要把面元離散密度控制在合適的范圍內(nèi),弦向50×展向40的離散密度能夠保證計算準(zhǔn)確。

圖8 不同面元離散密度下的推力

2.2 氣動性能

2.2.1 推力對比

選擇軸向來流風(fēng)速為15 m/s,轉(zhuǎn)速為4 000,4 500,5 000,5 500,6 000 r/min的工況進(jìn)行仿真,可獲得面元-渦粒子法和有限體積法推力。

如圖9所示,使用面元-渦粒子法與有限體積方法獲得的推力趨勢一致;同時在4 000 r/min時誤差為8%,在4 500 r/min以上時推力誤差降低。

圖9 不同轉(zhuǎn)速下螺旋槳推力大小

2.2.2 葉表壓力分布對比

選擇轉(zhuǎn)速為5 000 r/min,軸向前飛速度15 m/s的工況。在25%,50%,75%葉高處,對比面元-渦粒子法和有限體積法解出的壓力分布。

結(jié)合圖10,對比不同葉高處的靜壓分布發(fā)現(xiàn):

1) 在前緣滯止點處和吸、壓力面大部分區(qū)域內(nèi),2種方法能獲得基本一致的靜壓。

2) 在距離壓力面前緣0.1倍弦長位置處,使用面元-渦粒子法算得靜壓偏低。在無黏假設(shè)下,面元法求出的氣體加速作用更明顯,而真實氣體存在黏性,在葉背處加速作用稍弱,產(chǎn)生靜壓計算誤差。

圖10 25%葉高處螺旋槳壓力分布對比

2.2.3 尾跡流場

從圖11a)中可以看出,螺旋槳尾跡區(qū)域有明顯的周期性葉尖渦存在。趙寅宇[19]證實在葉尖渦粒子對誘導(dǎo)速度計算有主要作用,而在圖11a)中也能看到葉尖渦強(qiáng)度比槳盤其他展向位置的渦量要更高。圖11b)是在尾跡區(qū)內(nèi)周向速度大小,從圖中可以看出,每經(jīng)過高渦量區(qū)域會產(chǎn)生軸向加速作用。

圖11 垂直方向尾跡區(qū)的渦量和速度分布

為了進(jìn)一步揭示流場信息,下面對z軸坐標(biāo)為-0.4 m處的切面上渦量和誘導(dǎo)速度進(jìn)行分析。

為了提高數(shù)值計算收斂性,在使用有限體積法計算時考慮槳轂的影響。從圖12可以看出,螺旋槳葉尖半徑區(qū)域的尾跡區(qū)渦量仍然較高,呈現(xiàn)對稱分布,與面元-渦粒子法中的1-1旋渦對應(yīng),在槳轂半徑附近區(qū)域形成了較小的旋渦,也即是3-3旋渦。

圖12 出口水平切面上的尾跡渦量

為了進(jìn)一步驗證2種方法解出的尾跡速度分布一致性,提取對角線(圖12中黑色虛線)上速度分布，如圖13所示。

從圖13中可以看出,面元-渦粒子法和有限體積法解出的速度分布均呈現(xiàn)出“M”形,使用上述方法能夠獲得一致的尾跡速度分布。

圖13 距出口軸向0.4 m處橫截面軸向速度分布

2.3 聲學(xué)性能對比

在本節(jié)中將結(jié)合各倍頻下的聲壓級,對比面元-渦粒子法和有限體積法分別與Lowson聲學(xué)模型結(jié)合求出的聲學(xué)結(jié)果,分析螺旋槳聲學(xué)特性以及不同方法之間的差異。取垂直于螺旋槳槳盤平面布置觀測點,觀測點半徑為螺旋槳半徑的5倍。

圖14 不同葉片通過頻率下的總聲壓級

結(jié)合圖14可以看出,使用2種不同的方法解出的總聲壓級分布形式具有如下特征:

1) 使用2種方法解出的指向性一致,均呈現(xiàn)出以螺旋槳槳盤面為對稱面,“8”形的分布形式,具有典型的偶極子聲源分布特征[20];

2) 隨著倍頻的不斷增加,2種方法預(yù)測出的聲壓級均有明顯降低,但聲壓級指向性特征未發(fā)生改變,仍然軸向聲壓級最高,而在槳盤平面上最低。

上述現(xiàn)象說明,面元-渦粒子法能夠較好捕捉到聲壓級指向性特征,有效預(yù)測螺旋槳噪聲。

2.4 計算效率對比

為了評估面元-渦粒子法的計算效率,對比分別使用有限體積法與面元-渦粒子法計算相同時間步的耗時情況。數(shù)值計算在同一臺服務(wù)器上完成,服務(wù)器的CPU型號是英特爾銳龍Gold 6132,該CPU運(yùn)算主頻是2.6 GHz,運(yùn)算性能參數(shù)對比如表4所示。

表4 計算耗時對比

對比單核運(yùn)算耗時可以發(fā)現(xiàn),面元-渦粒子法計算效率明顯高于有限體積法。面元-渦粒子法能有效提高螺旋槳的非定常氣動和聲學(xué)仿真效率,節(jié)省計算資源,在設(shè)計階段具有較高的應(yīng)用價值。

3 結(jié) 論

本文發(fā)展了基于面元-渦粒子法的葉輪機(jī)械非定常氣動及噪聲的快速預(yù)測方法。使用該方法與二階有限體積法預(yù)測某型螺旋槳非定常氣動性能,發(fā)現(xiàn)在推力、葉表壓力以及尾跡速度分布等氣動方面和聲學(xué)方面具有較好的預(yù)測精度,研究結(jié)論如下:

1) 對比不同面元離散密度的計算結(jié)果,發(fā)現(xiàn)需要將面元尺寸控制在合適的范圍內(nèi),使之同時滿足能準(zhǔn)確描述葉表幾何尺寸和誘導(dǎo)作用準(zhǔn)確的需求。使用面元-渦粒子法可以準(zhǔn)確獲得螺旋槳的推力、葉表壓力以及尾跡速度分布,在不同葉高截面處壓力分布最大相對誤差為5%以內(nèi);

2) 渦粒子法在處理尾跡渦強(qiáng)發(fā)展時具有低的數(shù)值耗散優(yōu)勢。與二階有限體積法相比,使用面元-渦粒子法可以獲得明顯的葉尖渦誘導(dǎo)作用,并且由于不存在轉(zhuǎn)靜交界面的數(shù)值耗散,尾跡渦強(qiáng)連續(xù)性也隨之提高;

3) 使用面元-渦粒子法可以獲得與有限體積法指向性一致的聲壓級分布,1BPF下前向輻射角60°的聲壓級誤差為5%,該特征角度下,2BPF下誤差為8%。除此之外,面元-渦粒子法的計算效率更高,將該方法與聲學(xué)模型混合適用于螺旋槳非定常氣動噪聲性能的快速評估和進(jìn)一步優(yōu)化設(shè)計。