風(fēng)力機非定常氣動力降階模型

郭 強，劉鵬寅，竺曉程

（1.上海汽車集團股份有限公司乘用車公司，上海201804；2.上海交通大學(xué)機械與動力工程學(xué)院，上海200240）

眾所周知，準(zhǔn)確的非定常氣動力模型是風(fēng)力機氣彈設(shè)計的重要組成部分.而氣彈的不穩(wěn)定和顫振是風(fēng)力機設(shè)計必須考慮的問題，因此準(zhǔn)確預(yù)測風(fēng)力機葉片所受動態(tài)氣動載荷在風(fēng)力機設(shè)計中尤為重要.由于整個設(shè)計流程是一個不斷迭代的過程，因此所需的氣彈計算模型不僅要準(zhǔn)確而且要快速，尤其是對于氣彈現(xiàn)象復(fù)雜的風(fēng)力機葉片［1］.風(fēng)力機葉片的非定常氣動載荷主要是由于在自然環(huán)境中，大氣紊流、塔架影響、地面邊界層效應(yīng)以及偏航工況等因素影響，使得風(fēng)力機葉片處于非定常流動區(qū)域.在此情況下，葉素的有效攻角在不斷變化，并會導(dǎo)致動態(tài)氣動載荷的出現(xiàn).

這一動態(tài)氣動載荷的特征取決于邊界層是附著還是分離.如果葉片表面流動是附著的，翼型所受的動態(tài)氣動載荷與靜態(tài)氣動載荷差別不大，主要表現(xiàn)為一種延時現(xiàn)象.如果葉片表面流動出現(xiàn)分離，葉片振動將使得有效攻角呈非定常變化，同時導(dǎo)致動態(tài)失速的發(fā)生.在流動分離工況下，翼型的動態(tài)氣動載荷和靜態(tài)氣動載荷會有相當(dāng)大的區(qū)別，動態(tài)失速使葉片所受的最大氣動載荷大大增加.而隨著風(fēng)力機直徑的增大，過大的氣動載荷增加了風(fēng)力機穩(wěn)定控制的難度.從穩(wěn)定方面考慮，至少在升力預(yù)測時需要加入動態(tài)失速模型以減小預(yù)測誤差［2］.

從非定常氣動力預(yù)測精度來看，隨著計算機技術(shù)的顯著進步以及數(shù)值計算和湍流模型方面所取得的進展，通過CFD 方法來研究非定常氣動力的動態(tài)響應(yīng)變得越來越切實可行.大量這方面的工作［3-5］不僅對了解非定常氣動力的變化發(fā)展有積極作用，同時也證明了CFD 方法在常規(guī)風(fēng)力機和直升機槳葉設(shè)計中必會起到越來越重要的作用.但是較高精度的CFD 方法需要較高的計算成本.

因此，依靠少量數(shù)值計算算例結(jié)果建立預(yù)測模型，并對其他工況進行預(yù)測模擬的方法逐漸引起研究人員的興趣.降階模型方法是這類非定常氣動力預(yù)測模擬方法中較常用的一種，常用的2種降階模型有：ARMA 模型［6］和Volterra級數(shù)模型［7］.這類方法的主要優(yōu)點是其精度可隨著CFD 方法的改進而提升.當(dāng)流動復(fù)雜程度增加時，基于高適用性的CFD 方法的降階模型將比半經(jīng)驗?zāi)Ｐ透鼮榭煽?

由于三維問題的復(fù)雜程度較高，筆者將風(fēng)力機葉片的動態(tài)氣動載荷問題簡化為二維翼型進行振蕩運動時的氣動力變化.在附著流動和分離流動情況下，分別對翼型的升力系數(shù)采用2種不同降階模型進行預(yù)測，并將預(yù)測結(jié)果與CFD 計算結(jié)果進行對比，以驗證降階模型的有效性.

1 CFD計算

采用帶轉(zhuǎn)捩修正的k-ωSST 湍流模型.圖1為計算所用的網(wǎng)格和翼型表面附近的網(wǎng)格.采用三角形和四邊形混合的網(wǎng)格，整個計算區(qū)域分為外部靜止區(qū)域（簡稱外部區(qū)域）和內(nèi)部圍繞翼型振蕩的運動區(qū)域（簡稱內(nèi)部區(qū)域）.外部區(qū)域全部采用四邊形網(wǎng)格，內(nèi)部區(qū)域采用三角形和四邊形混合網(wǎng)格.為了減小界面通量傳遞的插值誤差，內(nèi)、外區(qū)域交界面具有相同大小和數(shù)量的一致網(wǎng)格.翼型周圍采用C 型四邊形邊界層網(wǎng)格，翼型表面網(wǎng)格節(jié)點總數(shù)為400.第一層網(wǎng)格的布置確保y＋在1.0量級處，且保證網(wǎng)格的法向膨脹率不大于1.1.外部進、出口邊界位于距翼型約20倍弦長處.入口邊界設(shè)置速度分量及合適的入流湍流度，出口邊界為靜壓力邊界，翼型表面的邊界條件為振蕩運動的無滑移壁面.經(jīng)過時間步長無關(guān)性檢查，采用的時間步長為0.001s，每個工況都至少經(jīng)過5個連續(xù)周期的計算以得到周期性穩(wěn)定的流場［4］.需要指出的是這樣一個非定常工況的CFD 計算需要十幾個小時，對于需要大量工況數(shù)據(jù)的設(shè)計和優(yōu)化工作而言計算成本較高.

圖1 計算網(wǎng)格Fig.1 Computational mesh

2 降階模型

2.1 ARMA 模型

ARMA 模型將某時刻動態(tài)氣動載荷的響應(yīng)認(rèn)為是當(dāng)前時刻的動態(tài)輸入和前幾個時刻的動態(tài)輸入與響應(yīng)的線性組合.因此當(dāng)前時刻的動態(tài)響應(yīng)可以表示為：

式中：u（n）為動態(tài)輸入，此處認(rèn)為是n時刻的攻角；y（n）為非定常狀態(tài)的動態(tài)氣動載荷響應(yīng)，在本文中響應(yīng)y為升力系數(shù)；na和nb分別為可自定義的模型的階數(shù)；ani和bni為待辨識的模型參數(shù).

將式（1）以矩陣形式給出，可得到ARMA 模型的另一表達(dá)式：

式中：n＝ns＋1，ns＋2，…，N；ns＝max［na，nb］；N為離散時間步的總數(shù).

采用最小二乘法進行辨識，可以得到模型的未定參數(shù)：

通過以上方法提取一個工況算例的數(shù)據(jù)作為樣本，通過辨識得到ARMA 模型的待定參數(shù)，從而建立ARMA 模型，并依靠建立的ARMA 模型對其他工況的動態(tài)響應(yīng)進行預(yù)測.

2.2 Volterra級數(shù)模型

Volterra級數(shù)模型是較常用的非線性響應(yīng)模型.對一個離散系統(tǒng)來說，輸入響應(yīng)可由式（7）表示：

式中：h0為穩(wěn)態(tài)響應(yīng)值；hi（n-k1，n-k2，…，n-ki）為第i階Volterra核.

由于Volterra核隨階數(shù)的變化呈指數(shù)增長，求解高階的Volterra核需要較高的計算資源確定，所以對于一個非線性系統(tǒng)，二階或三階Volterra級數(shù)是最常用的模型.筆者選用二階Volterra級數(shù)模型.

設(shè)定一個記憶長度M，二階的Volterra級數(shù)可以變化為

若將輸入以矩陣形式表示：

由此可以將二階的Volterra級數(shù)模型改寫為矩陣形式：

其中

這樣待辨識的模型參數(shù)就可由最小二乘法得到：

一旦辨識得到各階Volterra核，就可以建立Volterra級數(shù)模型并預(yù)測非定常氣動載荷.需要特別指出的是，以上2個模型的辨識所需時間在1～2 min內(nèi)，而模型建立后預(yù)測某一工況結(jié)果所需的時間在幾秒鐘左右.

3 不同工況的降階模型預(yù)測

采用S809翼型進行計算，其靜態(tài)失速攻角在10°左右.該翼型的升力系數(shù)在失速攻角（最高值）附近保持一段較高的數(shù)值，而阻力系數(shù)在很大的攻角變化范圍內(nèi)較小.這些特性使得翼型在變工況下保持較好的氣動特性.為了得到動態(tài)失速下的翼型非定常氣動載荷，定義風(fēng)力機翼型繞轉(zhuǎn)軸（x／c＝0.25）周期振蕩的變化形式為：

式中：α0為振蕩的平均攻角；Δα為振蕩的振幅；f為振蕩頻率.

同時定義衰減頻率為

式中：c為弦長；U為來流速度.

衰減頻率反映了沿翼型弦向的主流運動和繞轉(zhuǎn)軸的振蕩運動之間的相互作用，其值表征振蕩運動對主流運動影響的大小，所以當(dāng)給出振蕩的平均攻角、振幅和衰減頻率后，某一特定翼型的振蕩變化就可以確定.

所有模型的辨識數(shù)據(jù)都是由CFD 數(shù)值計算得到的.對每個非定常氣動計算結(jié)果，定義時間t0和tf為觀測的初始和終止時刻.對一個周期內(nèi)的數(shù)據(jù)進行處理后得到：

3.1 流動附著工況

將α0＝4°、k＝0.050和Δα＝2°工況下的CFD計算作為樣本，根據(jù)式（6）建立ARMA 模型.由于更高的階數(shù)已不能提高模型的精度，通過對比選取na＝2，nb＝1.當(dāng)翼型在這一小攻角（遠(yuǎn)小于靜態(tài)失速攻角）處振蕩時，流動未發(fā)生分離，因此升力系數(shù)的變化十分平緩，非線性現(xiàn)象很輕微.這種情況下采用模擬弱非線性系統(tǒng)的ARMA 模型進行動態(tài)氣動載荷的預(yù)測較有效.

圖2給出了樣本工況的升力系數(shù)隨攻角變化的氣動力遲滯回線.由圖2可知，與靜態(tài)不同，在附著流動時，同一攻角下升力系數(shù)在上仰和下俯過程中是不相同的.升力系數(shù)在整個翼型振蕩周期內(nèi)形成一個橢圓形的氣動力遲滯閉環(huán).由于流動處于附著狀態(tài)，所以動態(tài)升力系數(shù)與靜態(tài)升力系數(shù)在數(shù)值上差別不大.在非定常和定常升力之間存在一個時間延遲，這一時間延遲可以由Theodorsen理論計算得到.比較ARMA 模型、CFD 模擬和Theodorsen理論計算的結(jié)果可以看出，ARMA 模型對樣本工況的氣動力遲滯回線與CFD 計算結(jié)果基本吻合.同時Theodorsen理論計算結(jié)果與兩者的差別很小，說明CFD 計算結(jié)果與理論計算結(jié)果吻合良好.由于Theodorsen理論計算是在氣體不可壓縮的假設(shè)下建立的，而CFD 計算需要大量計算成本，所以采用ARMA 模型在計算成本較低的前提下獲得滿意的動態(tài)氣動載荷預(yù)測精度就具有普適性的優(yōu)勢.

圖2 不同模型的升力系數(shù)預(yù)測Fig.2 Lift coefficient predicted by different models

采用所建ARMA 模型對其他工況的非定常氣動載荷進行預(yù)測，注意攻角的變化范圍應(yīng)小于靜態(tài)失速攻角.圖3給出了ARMA 模型對平均攻角α0＝4°時，振蕩振幅Δα＝2°、1°和衰減頻率k＝0.050、0.060的4個工況以及平均攻角α0＝2°，Δα＝2°時衰減頻率k＝0.050、0.060的2個工況下升力系數(shù)的預(yù)測結(jié)果.從圖3 可以看出，對比CFD 計算結(jié)果，ARMA 模型的預(yù)測結(jié)果十分精確.ARMA 模型在附著流動工況下對非定常氣動載荷有很好的預(yù)測能力.

3.2 流動分離工況

當(dāng)攻角在較大數(shù)值（接近或超過靜態(tài)失速攻角）處變化時，由于分離的發(fā)生，非線性現(xiàn)象得到加強.圖4給出了CFD 計算得出的升力系數(shù)在α0＝8°、k＝0.050和Δα＝5.5°時的變化曲線，同時建立Volterra級數(shù)模型和ARMA 模型樣本工況進行對比分析.由于過高的記憶長度會提高模型系統(tǒng)辨識的難度，此處選取記憶長度M＝5.從圖4 可以看出，在流動分離工況下，升力系數(shù)遲滯回線不再呈橢圓形，這與流動附著時的Theodorsen理論計算結(jié)果不同.在翼型上仰階段，失速有被延遲的跡象，而在下俯階段重新附著的過程也存在遲滯.同時，升力系數(shù)的遲滯回線出現(xiàn)交叉，說明流動的復(fù)雜程度較高.而遲滯回線的閉環(huán)明顯要大于附著流動時，說明在流動分離工況下，動態(tài)氣動載荷的變化范圍大于附著流動工況.

圖3 ARMA 模型對不同工況下升力系數(shù)的預(yù)測Fig.3 Lift coefficient predicted by ARMA model under different conditions

圖4 ARMA 模型和Volterra級數(shù)模型的預(yù)測對比Fig.4 Comparison of prediction results between ARMA and Volterra series model

在流動分離工況下，ARMA 模型對CFD 計算結(jié)果的預(yù)測精度大大降低，表明此工況下ARMA模型很難捕捉到響應(yīng)的實時變化.這是因為ARMA模型只適用于模擬弱非線性的問題，而升力系數(shù)在流動分離時的動態(tài)響應(yīng)呈現(xiàn)較高的非線性特征.Volterra級數(shù)模型對較強非線性問題的模擬能力高于ARMA 模型，所以從圖4可以看出，Volterra級數(shù)模型對于訓(xùn)練樣本工況的預(yù)測與CFD 計算結(jié)果吻合更好.采用Volterra級數(shù)模型對其他工況進行預(yù)測，其結(jié)果與CFD 計算結(jié)果的對比見圖5.圖中所示的測試工況分別是Δα＝5.5°、10°和k＝0.026，0.050，0.060，0.077的4種組合.由圖5可知，Volterra級數(shù)模型對樣本工況以及其他工況除細(xì)微差別外整體預(yù)測精度很高.說明Volterra級數(shù)模型在流動分離工況下可以比較精確地預(yù)測翼型所受的非定常氣動載荷.圖4和圖5的計算結(jié)果也表明，記憶長度M＝5已經(jīng)能夠保證模型的精度.

圖5 Volterra級數(shù)模型對不同工況下升力系數(shù)的預(yù)測Fig.5 Lift coefficient predicted by Volterra series model under different conditions

4 結(jié) 論

（1）ARMA 模型在流動附著工況下可以很好地預(yù)測翼型非定常氣動載荷的變化，但在流動分離工況下其預(yù)測能力較差.

（2）Volterra級數(shù)模型在流動分離工況下的預(yù)測能力較ARMA 模型有了很大提高.其對非定常氣動載荷的預(yù)測結(jié)果與CFD 計算結(jié)果吻合較好，相比具有較高的精度.

（3）采用ARMA 模型和Volterra級數(shù)模型預(yù)測各自適應(yīng)工況的非定常氣動載荷是可行的.由于ARMA 模型和Volterra級數(shù)模型只需要一個工況的CFD 計算數(shù)據(jù)作為樣本去辨識各自的模型待定參數(shù)，而一旦辨識好后只需幾秒時間就可預(yù)測其他工況的非定常氣動載荷.因此，這2種只需較小計算成本的降階模型為非定常氣動載荷的預(yù)測、相關(guān)的氣動彈性分析以及風(fēng)力機葉片的設(shè)計和優(yōu)化提供了一種有效的方法.

［1］LEISHMAN J G，NGUYEN K Q.State-space representation of unsteady airfoil behavior［J］.AIAA Journal，1990，28（5）：836-844.

［2］HANSEN MOL.Aerodynamics of wind turbines［M］.London，UK：Earthscan，2008.

［3］劉雄，馬新穩(wěn)，沈世，等.風(fēng)力機柔性葉片振動變形對其氣動阻尼的影響分析［J］.空氣動力學(xué)學(xué)報，2013，31（3）：407-412.

LIU Xiong，MA Xinwen，SHEN Shi，etal.Analysis of the influence of vibration and deformation of the blade on the aerodynamic damping［J］.Acta Aerodynamica Sinica，2013，31（3）：407-412.

［4］俞國華.水平軸風(fēng)力機葉片失速問題研究［D］.上海：上海交通大學(xué)，2013.

［5］周正，李春.風(fēng)力機翼型等速上仰動態(tài)失速數(shù)值模擬［J］.能源研究與信息，2013，29（4）：196-200.

ZHOU Zheng，LI Chun.Numerical simulation of identical pitching wind turbine airfoil dynamic stall［J］.Energy Research and Information，2013，29（4）：196-200.

［6］楊錫運，孫翰墨.基于時間序列模型的風(fēng)電場風(fēng)速預(yù)測研究［J］.動力工程學(xué)報，2011，31（3）：203-208.

YANG Xiyun，SUN Hanmo.Wind speed prediction in wind farms based on time series model［J］.Journal of Chinese Society of Power Engineering，2011，31（3）：203-208.

［7］徐敏，陳剛，陳士櫓，等.基于非定常氣動力低階模型的氣動彈性主動控制律設(shè)計［J］.西北工業(yè)大學(xué)學(xué)報，2005，22（6）：748-752.

XU Min，CHEN Gang，CHEN Shilu，etal.Design method of active control law based on reduced order model of unsteady aerodynamics for aeroelasticity［J］.Journal of Northwestern Polytechnical University，2005，22（6）：748-752.