風(fēng)力機(jī)翼型非定常流場(chǎng)POD和EPOD分析

孫 翀， 田 甜， 竺曉程， 杜朝輝

(上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院， 上海 200240)

風(fēng)能是當(dāng)前最具大規(guī)模開(kāi)發(fā)價(jià)值的可再生能源之一，風(fēng)力機(jī)作為風(fēng)能利用的主要裝置，在向大功率發(fā)展的同時(shí)，其產(chǎn)生的噪聲對(duì)人們生活的影響也日益明顯，降噪已經(jīng)成為風(fēng)力機(jī)研究的重要方面[1-2].二維翼型作為風(fēng)力機(jī)葉片最基本單元，其非定常流場(chǎng)對(duì)風(fēng)力機(jī)整機(jī)性能以及氣動(dòng)噪聲水平具有重要的影響[3].在氣動(dòng)噪聲方面，有研究表明由翼型繞流產(chǎn)生的氣動(dòng)噪聲包括：尾緣噪聲、鈍尾噪聲和失速噪聲[4].其中，現(xiàn)代大型水平軸風(fēng)力機(jī)氣動(dòng)噪聲又以渦團(tuán)通過(guò)葉片尾緣后形成的尾緣噪聲為主[5].為此，需對(duì)翼型的非定常流動(dòng)進(jìn)行精細(xì)的求解和分析.

實(shí)驗(yàn)測(cè)試和計(jì)算流體力學(xué)(CFD)是研究翼型非定常流動(dòng)的主要方法.其中，CFD方法擁有花費(fèi)低，限制少等優(yōu)勢(shì)，因此被廣泛運(yùn)用于翼型流場(chǎng)數(shù)據(jù)的獲取.但是，由CFD得到的非定常流場(chǎng)數(shù)據(jù)非常龐大，直接對(duì)其主要流動(dòng)結(jié)構(gòu)分析較為困難.降階模型方法通過(guò)一組低維變量構(gòu)成的特征模態(tài)來(lái)表示非定常流場(chǎng).其中，本征正交分解(POD)方法是模態(tài)分析方法中較為傳統(tǒng)和常用的一種方法，通過(guò)POD方法可以提取獲得非定常流場(chǎng)的主要流動(dòng)模態(tài)，其特征值直接表征對(duì)應(yīng)模態(tài)的能量[6].在翼型非定常流動(dòng)的POD應(yīng)用方面，董圣華等[7]研究了超臨界翼型在跨聲速抖振下的非定常流場(chǎng)，采用POD方法提取了導(dǎo)致抖振的主要流動(dòng)模態(tài).Zhao等[8]使用POD方法研究了凹凸前緣翼型流動(dòng)控制機(jī)理和氣動(dòng)特性.

對(duì)于流場(chǎng)中某個(gè)特定的區(qū)域或者某些觀測(cè)變量的研究，若其能量占流動(dòng)總能量的比例較小時(shí)，POD方法會(huì)有一定的局限性， Borée[9]提出了擴(kuò)展本征正交分解(EPOD)方法，可以獲得與某個(gè)區(qū)域中與觀測(cè)變量最相關(guān)的流動(dòng)結(jié)構(gòu).Schlegel等[10]基于EPOD方法，研究了射流中于氣動(dòng)噪聲最相關(guān)的流動(dòng)結(jié)構(gòu).

本文通過(guò)大渦模擬(LES)方法，數(shù)值計(jì)算風(fēng)力機(jī)S809翼型在不同攻角下的非定常流場(chǎng)，并求解Ffowcs Williams-Hawkings(FW-H)方程獲得遠(yuǎn)場(chǎng)翼型遠(yuǎn)場(chǎng)氣動(dòng)噪聲.通過(guò)POD方法對(duì)流場(chǎng)數(shù)據(jù)進(jìn)行降階處理，獲得主要的非定常流動(dòng)結(jié)構(gòu)，并進(jìn)一步通過(guò)EPOD方法，構(gòu)建遠(yuǎn)場(chǎng)氣動(dòng)噪聲與翼型非定常流動(dòng)的關(guān)系，提取與氣動(dòng)噪聲最相關(guān)的流動(dòng)結(jié)構(gòu)，揭示風(fēng)力機(jī)翼型氣動(dòng)噪聲形成和傳播的物理機(jī)制，為大型風(fēng)力機(jī)降噪提供理論依據(jù).

1 數(shù)值計(jì)算

對(duì)典型風(fēng)力機(jī)S809翼型進(jìn)行數(shù)值計(jì)算，翼型弦長(zhǎng)c=0.2 m，計(jì)算域原點(diǎn)為翼型氣動(dòng)中心(1/4弦長(zhǎng)位置)，半徑為20c，展向長(zhǎng)度為c.使用ICEM (The Integrated Computer Engineering and Manufacturing)軟件對(duì)計(jì)算域做結(jié)構(gòu)化網(wǎng)格劃分，周向288個(gè)節(jié)點(diǎn)，徑向215個(gè)節(jié)點(diǎn)，展向41個(gè)節(jié)點(diǎn)，總網(wǎng)格數(shù)達(dá)到250萬(wàn)，全局網(wǎng)格如圖1所示.為提升壁面邊界層內(nèi)的計(jì)算精度，對(duì)翼型附近網(wǎng)格進(jìn)行加密，近壁面第1層網(wǎng)格高度設(shè)置為0.005 mm，該值保證了壁面無(wú)量綱距離y+小于1，并將網(wǎng)格壁法向膨脹比設(shè)置為1.1以實(shí)現(xiàn)邊界層計(jì)算高分辨率.

計(jì)算邊界條件參考文獻(xiàn)[11-12]，如圖2所示.進(jìn)出口設(shè)置于翼型前后20c處，進(jìn)口給定來(lái)流速度v0=38.63 m/s，出口背壓設(shè)置為一個(gè)大氣壓，由于風(fēng)力機(jī)翼型來(lái)流馬赫數(shù)低，計(jì)算域足夠大，該邊界設(shè)置與遠(yuǎn)場(chǎng)條件非常接近；展向邊界設(shè)置為對(duì)稱(chēng)面；翼型表面設(shè)置為無(wú)滑移壁面.使用商業(yè)軟件ANSYS CFX進(jìn)行數(shù)值計(jì)算，湍流模擬使用LES Wale模型.計(jì)算時(shí)間步長(zhǎng)設(shè)置為dτ=0.02 ms，使計(jì)算平均庫(kù)朗數(shù)在1左右，以保證對(duì)湍流發(fā)展的數(shù)值模擬精度.

圖1 翼型計(jì)算網(wǎng)格Fig.1 Computation mesh of airfoil

圖2 計(jì)算域邊界Fig.2 Computation domain

2 POD及EPOD方法

(1)

j=1,2,…,n

式中：φj(ξ)為第j個(gè)空間模態(tài)；aj(t)為第j個(gè)時(shí)間系數(shù).POD方法則是為這樣的分解形式尋找一組正交的模態(tài)[13].為尋找這樣的一組正交模態(tài)，可以采用快照方法.首先，定義離散時(shí)刻下的流場(chǎng)快照zi∈Rm×1，m為流動(dòng)變量q的空間維數(shù)，則有：

(2)

i=1,2,…,n

式中：ti為研究時(shí)間t內(nèi)的離散時(shí)刻.將這些快照組成快照矩陣：

(3)

求解快照矩陣協(xié)矩陣的特征值問(wèn)題：

ZTZΨj=λjΨj

(4)

j=1,2,…,n

式中：ψj為第j個(gè)特征向量；λj為第j個(gè)特征值.該特征值λj即為流動(dòng)變量q的POD特征值λqj，而q的POD特征模態(tài)φPOD,j由特征向量ψj的變換獲得：

(5)

j=1,2,…,n

動(dòng)能可以通過(guò)POD特征值之和來(lái)表征，每個(gè)模態(tài)能量占總能量的比重可以由對(duì)應(yīng)特征值與特征值之和的比值表示.因此，將特征值從大到小排列，所對(duì)應(yīng)的前幾階POD模態(tài)便為表征流場(chǎng)波動(dòng)的主要模態(tài).

根據(jù)式(1)，通過(guò)POD特征模態(tài)φPOD,j的原始快照zi為

(6)

i=1,2,…,n

式中：βqi,j為模態(tài)系數(shù).

通過(guò)求解式(6)可獲得模態(tài)系數(shù)矩陣：

(7)

(8)

若有另一個(gè)時(shí)空域Ω上的流動(dòng)參數(shù)l，則EPOD研究了l對(duì)q的模態(tài)投影的關(guān)系.為獲得EPOD模態(tài)，首先對(duì)參數(shù)q做POD分解，可以得到參數(shù)l的POD特征值λl，以及參數(shù)l的POD模態(tài)系數(shù)βl,則q關(guān)于l的EPOD模態(tài)φEP滿足如下關(guān)系[9]：

(9)

EPOD模態(tài)被證明為表示q中所有與l相關(guān)的流動(dòng)結(jié)構(gòu)[9]，由模態(tài)特征值λl的大小排序，可以得到與l最相關(guān)的流動(dòng)結(jié)構(gòu).

3 計(jì)算結(jié)果

圖3 翼型表面壓力系數(shù)Fig.3 Pressure coefficients of airfoil surface

根據(jù)以上數(shù)值計(jì)算結(jié)果，通過(guò)求解FW-H方程，可以獲得聲壓監(jiān)測(cè)點(diǎn)下的翼型氣動(dòng)噪聲.當(dāng)翼型攻角為8° 時(shí)，氣動(dòng)中心正上方2 m位置處的監(jiān)測(cè)點(diǎn)聲壓頻譜如圖4(a)所示.其中：f為頻率；SPL為聲壓級(jí).氣動(dòng)噪聲呈寬頻凸起特性，駝峰對(duì)應(yīng)頻率在500～800 Hz范圍，整體變化規(guī)律與NREL的BPM(Brooks, Pope, Marcolin)翼型自噪聲預(yù)測(cè)模型結(jié)果接近.

圖4 翼型氣動(dòng)噪聲Fig.4 Aerodynamic noise of airfoil

為研究翼型氣動(dòng)噪聲指向性，需要計(jì)算不同傳播方向上的氣動(dòng)噪聲強(qiáng)度，翼型氣動(dòng)中心為圓點(diǎn)，以10c(即2 m)為半徑，圓心角每隔5° 布置一個(gè)聲壓監(jiān)測(cè)點(diǎn)，共72個(gè)監(jiān)測(cè)點(diǎn)。計(jì)算這些監(jiān)測(cè)點(diǎn)上的時(shí)均總聲壓級(jí)，獲得的不同攻角下翼型氣動(dòng)噪聲聲壓級(jí)指向性分布如圖4(b)所示.由圖4(b)可知，當(dāng)翼型攻角為2° 和8° 時(shí)，氣動(dòng)噪聲強(qiáng)度在翼型上下對(duì)稱(chēng)，而在翼型前后明顯衰減，具有偶極子特征.翼型前部氣動(dòng)噪聲強(qiáng)度明顯弱于尾部，但隨著攻角的升高，翼型前部氣動(dòng)噪聲逐漸增強(qiáng).進(jìn)一步比較可以發(fā)現(xiàn)，翼型在攻角由8° 上升至14° 后，翼型吸力面?zhèn)仍肼曪@著增強(qiáng)，而壓力面?zhèn)仍肼晱?qiáng)度變化不大.

4 流場(chǎng)POD和EPOD分析

根據(jù)翼型在8°攻角下的非定常流場(chǎng)，每隔10步提取一個(gè)渦量流場(chǎng)數(shù)據(jù)作為一個(gè)快照樣本(快照間隔為2×10-4s )，共550個(gè)快照構(gòu)成快照矩陣Z，并做POD分解，獲得的POD特征值λ如圖5所示，其中：N為模態(tài)序號(hào).由圖5可知，POD特征值大小迅速下降，前4階POD模態(tài)占流場(chǎng)波動(dòng)總能量的45%，前40階模態(tài)占總能量的90%.

圖5 POD特征值Fig.5 Eigenvalues of POD

因此，POD方法可以有效地降階高維非定常流場(chǎng).分解得到的POD前4階模態(tài)及其對(duì)應(yīng)模態(tài)時(shí)間系數(shù)的頻譜分布如圖6所示.其中：x，y為翼型平面幾何坐標(biāo)；A為振幅.從圖6中可以看到，第1階模態(tài)結(jié)構(gòu)主要出現(xiàn)在尾緣位置.因此，翼型渦量流場(chǎng)最主要的非定常結(jié)構(gòu)表現(xiàn)為翼型附近尤其是尾緣處的渦團(tuán)，其模態(tài)系數(shù)的頻譜呈寬頻特征，且主要集中在中低頻(小于1 kHz).第2、3階模態(tài)結(jié)構(gòu)從前緣開(kāi)始，主要表現(xiàn)在翼型中部位置，該位置接近轉(zhuǎn)捩位置.在轉(zhuǎn)捩位置之前， 層流流動(dòng)附著壁面，隨著逆壓梯度逐漸增強(qiáng)，流動(dòng)開(kāi)始失穩(wěn)，層流邊界層脫離壁面，隨后湍流邊界層再次附著，在翼型壁面上形成層流分離泡，并在此處發(fā)生流動(dòng)轉(zhuǎn)捩.根據(jù)模態(tài)系數(shù)頻譜，該兩階模態(tài)有明顯的主要頻率 1.265 kHz.第4階POD模態(tài)類(lèi)似第1階模態(tài)，表現(xiàn)在轉(zhuǎn)捩位置之后的翼型湍流邊界層及尾緣附近的渦團(tuán)結(jié)構(gòu)，其尺度較大，頻率較低.

圖6 POD模態(tài)Fig.6 Modes of POD

為研究與氣動(dòng)噪聲形成并向下游傳播有關(guān)的流動(dòng)結(jié)構(gòu)，測(cè)點(diǎn)參考文獻(xiàn)[10]，在翼型正上方10c處，平行于來(lái)流方向，向下游每0.05 m布置一個(gè)聲壓監(jiān)測(cè)點(diǎn).提取這些監(jiān)測(cè)點(diǎn)的時(shí)域氣動(dòng)聲壓計(jì)算數(shù)據(jù)，并分析翼型非定常渦量場(chǎng)關(guān)于該列測(cè)點(diǎn)氣動(dòng)聲壓的EPOD，獲得的EPOD模態(tài)特征值如圖7所示.

由圖7可以看到，EPOD特征值大小迅速衰減，前4階EPOD模態(tài)就占聲壓波動(dòng)能量的75%，相對(duì)于氣動(dòng)噪聲，EPOD模態(tài)具有更好的降階特性.

前4階聲壓EPOD模態(tài)及其對(duì)應(yīng)的模態(tài)系數(shù)頻譜如圖8所示.將得到的各階EPOD模態(tài)單位向量化以與POD模態(tài)進(jìn)行對(duì)比.由圖8可以看到，聲壓EPOD模態(tài)不僅表現(xiàn)在翼型附近，更表現(xiàn)為尾跡中的渦團(tuán)結(jié)構(gòu)，可以發(fā)現(xiàn)經(jīng)過(guò)尾緣后衰減，再逐漸放大、飽和后再衰減的類(lèi)似波包的結(jié)構(gòu)，該結(jié)構(gòu)在射流噪聲中被認(rèn)為與氣動(dòng)噪聲源有關(guān)[15].通過(guò)EPOD模態(tài)系數(shù)頻譜可以看到，前幾階EPOD模態(tài)系數(shù)頻譜具有寬頻凸起的特征.第3、4階EPOD模態(tài)系數(shù)在頻率0.1～1.5 kHz范圍內(nèi)都具有較高的幅值，表明與氣動(dòng)噪聲生成和傳播的湍流結(jié)構(gòu)頻帶更寬，波動(dòng)更復(fù)雜.

圖7 EPOD特征值Fig.7 Eigenvalues of EPOD

進(jìn)一步分析翼型氣動(dòng)噪聲源，不同攻角下翼型聲壓EPOD模態(tài)如圖9所示.由圖9可知，當(dāng)翼型攻角為2° 時(shí)，聲壓EPOD模態(tài)在尾跡區(qū)的結(jié)構(gòu)與翼型攻角為8° 時(shí)的情況較接近，表現(xiàn)出翼型前部氣動(dòng)噪聲強(qiáng)度明顯較尾部弱的特征(見(jiàn)圖4(b))，說(shuō)明在小攻角下，翼型尾部和尾跡的渦量波動(dòng)主導(dǎo)氣動(dòng)噪聲.而在翼型攻角為14° 時(shí)，分離渦從前緣開(kāi)始脫落.通過(guò)EPOD模態(tài)可以發(fā)現(xiàn)，上方吸力面?zhèn)榷即嬖诓煌叨鹊摹⑴c氣動(dòng)噪聲相關(guān)的湍流結(jié)構(gòu)，并一直延伸至尾跡，而壓力面上模態(tài)結(jié)構(gòu)卻相對(duì)較弱，與翼型攻角為2° 和8° 時(shí)的EPOD模態(tài)結(jié)構(gòu)差異較大.在此攻角下，翼型上方吸力面?zhèn)鹊臍鈩?dòng)噪聲強(qiáng)度顯著高于下方.

圖8 EPOD模態(tài)Fig.8 Modes of EPOD

圖9 不同攻角下的EPOD模態(tài)Fig.9 EPOD modes at different angles of attack

5 結(jié)論

本文采用大渦模擬方法對(duì)風(fēng)力機(jī)S809翼型的非定常流場(chǎng)進(jìn)行了數(shù)值模擬，并結(jié)合FW-H方法獲得了翼型氣動(dòng)噪聲，采用POD提取了非定常渦量流場(chǎng)的主要降階模態(tài)結(jié)構(gòu)，并采用EPOD方法，揭示了與翼型氣動(dòng)噪聲最相關(guān)的非定常流動(dòng)模態(tài)結(jié)構(gòu)，獲得如下主要結(jié)論.

(1) 采用大渦模擬對(duì)翼型流動(dòng)的計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)試值基本吻合，計(jì)算得到的翼型氣動(dòng)噪聲在頻譜圖上呈寬頻凸起特征，峰值頻率在500～800 Hz范圍.在小攻角情況下，翼型氣動(dòng)噪聲指向性具有偶極子特征，隨著攻角的升高，翼型前方以及吸力面?zhèn)鹊臍鈩?dòng)噪聲強(qiáng)度有所增強(qiáng).

(2) POD方法可以有效地降階翼型非定常流場(chǎng)，前4階POD模態(tài)占流場(chǎng)波動(dòng)總能量的45%，模態(tài)結(jié)構(gòu)均出現(xiàn)在翼型表面的中后部分，表明非定常流動(dòng)的主要特征為翼型表面層流分離泡及尾緣附近渦團(tuán).

(3) 基于翼型氣動(dòng)噪聲的非定常流場(chǎng)EPOD分析，前4階EPOD模態(tài)占聲壓能量的75%，對(duì)氣動(dòng)噪聲降階較好.翼型尾跡中的湍流結(jié)構(gòu)都被發(fā)現(xiàn)與氣動(dòng)噪聲相關(guān)，表明尾緣噪聲是翼型最主要的氣動(dòng)噪聲源.而在大攻角下，吸力面分離中的湍流結(jié)構(gòu)也與氣動(dòng)噪聲相關(guān).