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

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

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

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

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

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

1 理論方法

1.1 面元法

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

2φ=0

(1)

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

(2)

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

(3)

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

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

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

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

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

圖1 NACA0012的余弦離散

1.2 渦粒子法

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

(7)

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

(8)

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

(9)

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

(10)

(11)

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

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

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

(16)

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

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

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

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

(17)

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

1.3 Lowson聲學(xué)模型

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

cn=an+ibn=

(18)

將時(shí)間微分項(xiàng)替換,引入周向和軸向力分量,分部積分后得到Lowson公式

(19)

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

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

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

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

(22)

1.4 面元-渦粒子法

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

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

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

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

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

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

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

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

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

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

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

3BPF對(duì)應(yīng)的聲波波長是0.243 m,為了能夠嚴(yán)格滿足點(diǎ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所示,為了驗(yàn)證該算例的準(zhǔn)確性,需要進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證:

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

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

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

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

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

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

2.1.2 面元離散尺度

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

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

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

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

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

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

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

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

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

2.2 氣動(dòng)性能

2.2.1 推力對(duì)比

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

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

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

2.2.2 葉表壓力分布對(duì)比

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

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

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

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

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

2.2.3 尾跡流場(chǎng)

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

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

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

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

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

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

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

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

2.3 聲學(xué)性能對(duì)比

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

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

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

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

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

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

2.4 計(jì)算效率對(duì)比

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

表4 計(jì)算耗時(shí)對(duì)比

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

3 結(jié) 論

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

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

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

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