風力機翼型非定常流場POD和EPOD分析

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

(上海交通大學 機械與動力工程學院， 上海 200240)

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

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

對于流場中某個特定的區域或者某些觀測變量的研究，若其能量占流動總能量的比例較小時，POD方法會有一定的局限性， Borée[9]提出了擴展本征正交分解(EPOD)方法，可以獲得與某個區域中與觀測變量最相關的流動結構.Schlegel等[10]基于EPOD方法，研究了射流中于氣動噪聲最相關的流動結構.

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

1 數值計算

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

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

圖1 翼型計算網格Fig.1 Computation mesh of airfoil

圖2 計算域邊界Fig.2 Computation domain

2 POD及EPOD方法

(1)

j=1,2,…,n

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

(2)

i=1,2,…,n

式中：ti為研究時間t內的離散時刻.將這些快照組成快照矩陣：

(3)

求解快照矩陣協矩陣的特征值問題：

ZTZΨj=λjΨj

(4)

j=1,2,…,n

式中：ψj為第j個特征向量；λj為第j個特征值.該特征值λj即為流動變量q的POD特征值λqj，而q的POD特征模態φPOD,j由特征向量ψj的變換獲得：

(5)

j=1,2,…,n

動能可以通過POD特征值之和來表征，每個模態能量占總能量的比重可以由對應特征值與特征值之和的比值表示.因此，將特征值從大到小排列，所對應的前幾階POD模態便為表征流場波動的主要模態.

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

(6)

i=1,2,…,n

式中：βqi,j為模態系數.

通過求解式(6)可獲得模態系數矩陣：

(7)

(8)

若有另一個時空域Ω上的流動參數l，則EPOD研究了l對q的模態投影的關系.為獲得EPOD模態，首先對參數q做POD分解，可以得到參數l的POD特征值λl，以及參數l的POD模態系數βl,則q關于l的EPOD模態φEP滿足如下關系[9]：

(9)

EPOD模態被證明為表示q中所有與l相關的流動結構[9]，由模態特征值λl的大小排序，可以得到與l最相關的流動結構.

3 計算結果

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

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

圖4 翼型氣動噪聲Fig.4 Aerodynamic noise of airfoil

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

4 流場POD和EPOD分析

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

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

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

圖6 POD模態Fig.6 Modes of POD

為研究與氣動噪聲形成并向下游傳播有關的流動結構，測點參考文獻[10]，在翼型正上方10c處，平行于來流方向，向下游每0.05 m布置一個聲壓監測點.提取這些監測點的時域氣動聲壓計算數據，并分析翼型非定常渦量場關于該列測點氣動聲壓的EPOD，獲得的EPOD模態特征值如圖7所示.

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

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

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

進一步分析翼型氣動噪聲源，不同攻角下翼型聲壓EPOD模態如圖9所示.由圖9可知，當翼型攻角為2° 時，聲壓EPOD模態在尾跡區的結構與翼型攻角為8° 時的情況較接近，表現出翼型前部氣動噪聲強度明顯較尾部弱的特征(見圖4(b))，說明在小攻角下，翼型尾部和尾跡的渦量波動主導氣動噪聲.而在翼型攻角為14° 時，分離渦從前緣開始脫落.通過EPOD模態可以發現，上方吸力面側都存在不同尺度的、與氣動噪聲相關的湍流結構，并一直延伸至尾跡，而壓力面上模態結構卻相對較弱，與翼型攻角為2° 和8° 時的EPOD模態結構差異較大.在此攻角下，翼型上方吸力面側的氣動噪聲強度顯著高于下方.

圖8 EPOD模態Fig.8 Modes of EPOD

圖9 不同攻角下的EPOD模態Fig.9 EPOD modes at different angles of attack

5 結論

本文采用大渦模擬方法對風力機S809翼型的非定常流場進行了數值模擬，并結合FW-H方法獲得了翼型氣動噪聲，采用POD提取了非定常渦量流場的主要降階模態結構，并采用EPOD方法，揭示了與翼型氣動噪聲最相關的非定常流動模態結構，獲得如下主要結論.

(1) 采用大渦模擬對翼型流動的計算結果與實驗測試值基本吻合，計算得到的翼型氣動噪聲在頻譜圖上呈寬頻凸起特征，峰值頻率在500～800 Hz范圍.在小攻角情況下，翼型氣動噪聲指向性具有偶極子特征，隨著攻角的升高，翼型前方以及吸力面側的氣動噪聲強度有所增強.

(2) POD方法可以有效地降階翼型非定常流場，前4階POD模態占流場波動總能量的45%，模態結構均出現在翼型表面的中后部分，表明非定常流動的主要特征為翼型表面層流分離泡及尾緣附近渦團.

(3) 基于翼型氣動噪聲的非定常流場EPOD分析，前4階EPOD模態占聲壓能量的75%，對氣動噪聲降階較好.翼型尾跡中的湍流結構都被發現與氣動噪聲相關，表明尾緣噪聲是翼型最主要的氣動噪聲源.而在大攻角下，吸力面分離中的湍流結構也與氣動噪聲相關.