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

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

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

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

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

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

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

1 理論方法

1.1 面元法

面元法以流體連續性為理論基礎,求解三維流場中勢函數:

2φ=0

(1)

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

(2)

根據固體壁面、尾跡面元不可穿透邊界條件,無限遠處勢函數誘導速度為0假設

(3)

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

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

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

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

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

圖1 NACA0012的余弦離散

1.2 渦粒子法

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

(7)

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

(8)

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

(9)

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

(10)

(11)

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

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

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

(16)

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

表1 常用的光滑函數和對應的Green函數

圖2 3種不同精度的Green函數

在表1常用的光滑函數和對應的Green函數中有

(17)

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

1.3 Lowson聲學模型

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

cn=an+ibn=

(18)

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

(19)

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

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

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

對第n階旋轉偶極子聲源的自由遠場聲壓表達式有

(22)

1.4 面元-渦粒子法

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

圖3 尾跡面元轉化為渦粒子

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

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

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

圖5 CFD和聲學模型耦合

2 數值計算結果

2.1 前處理參數設置

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

2.1.1 有限體積法網格尺度

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

3BPF對應的聲波波長是0.243 m,為了能夠嚴格滿足點聲源假設,網格尺度需要小于1/6聲波波長,所以網格尺度應該小于0.040 5 m。以表2中的網格尺度參數作為基準。

表2 網格大小設置參數 m

劃分多面體網格如圖6所示,為了驗證該算例的準確性,需要進行網格無關性驗證:

1) 邊界層網格高度

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

2) 尾跡區網格密度

尾跡區域中,在原網格基礎上疊加5%基礎尺寸網格,求得推力為11.474 N,與基準算例11.437 N相對誤差為0.32%,認為尾跡區網格密度滿足要求。

圖6 旋轉域網格示意圖

分別對邊界層、葉表、尾跡區進行網格加密并求解,證實使用該算例網格可以獲得準確結果。

2.1.2 面元離散尺度

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

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

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

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

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

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

3) 面元長寬比接近1時,計算結果更準確。

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

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

2.2 氣動性能

2.2.1 推力對比

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

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

圖9 不同轉速下螺旋槳推力大小

2.2.2 葉表壓力分布對比

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

結合圖10,對比不同葉高處的靜壓分布發現:

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

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

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

2.2.3 尾跡流場

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

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

為了進一步揭示流場信息,下面對z軸坐標為-0.4 m處的切面上渦量和誘導速度進行分析。

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

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

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

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

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

2.3 聲學性能對比

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

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

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

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

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

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

2.4 計算效率對比

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

表4 計算耗時對比

對比單核運算耗時可以發現,面元-渦粒子法計算效率明顯高于有限體積法。面元-渦粒子法能有效提高螺旋槳的非定常氣動和聲學仿真效率,節省計算資源,在設計階段具有較高的應用價值。

3 結 論

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

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

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

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