翼型擺角對(duì)氣動(dòng)性能的影響分析

季 康， 李 春，2， 葉 舟，2， 楊 陽(yáng)

（1.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093；2.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093）

風(fēng)力機(jī)在目前已建成的旋轉(zhuǎn)機(jī)械中直徑最大，然而隨著單機(jī)功率的增加，其葉片長(zhǎng)度也隨之越來越長(zhǎng)，這就會(huì)引起風(fēng)輪葉片的彎曲變形，最終導(dǎo)致氣動(dòng)性能的變化［1］.研究表明，當(dāng)風(fēng)力機(jī)葉片長(zhǎng)度增加至30m 以上時(shí)，葉片會(huì)產(chǎn)生柔性特性，再加上葉片本身材料的強(qiáng)度較高且模量較低，其幾何非線性特性將會(huì)明顯變強(qiáng)；而當(dāng)葉片長(zhǎng)度增加至60 m 以上時(shí)，尾緣由于同時(shí)受到氣動(dòng)力和結(jié)構(gòu)力，在時(shí)域范圍內(nèi)會(huì)出現(xiàn)截然不同的非定常特性［2－3］.

近年來，對(duì)翼型尾緣擺角的研究主要集中在飛行器方面，然而對(duì)于風(fēng)力機(jī)方面的研究相對(duì)偏向于葉片材料及控制等方面［4－9］.1999年，Andrew 等［10］第一次將主動(dòng)控制變形葉片技術(shù)應(yīng)用在風(fēng)力機(jī)葉片上，其對(duì)變形葉片的研究重點(diǎn)在于不同彎扭結(jié)合的葉片設(shè)計(jì)方式，而不改變?nèi)~片的整體結(jié)構(gòu)，文獻(xiàn)［11－14］給出了新型適應(yīng)性葉片，主要將葉片展向進(jìn)行了適當(dāng)?shù)膹澟ぴO(shè)計(jì)，但是并未涉及翼型尾緣變形后引發(fā)的氣動(dòng)效應(yīng).

本文選取NREL S809 翼 型［15］，通 過 數(shù) 值 模 擬商用軟件研究不同攻角、來流風(fēng)速和尾緣角度對(duì)其氣動(dòng)性能的影響，并對(duì)比變形翼型與原始翼型的升阻力系數(shù)及流線圖，分析了變形翼型上下表面壓力分布的變化規(guī)律.

1 尾緣變形翼型創(chuàng)建原理

對(duì)S809翼型尾緣進(jìn)行變形，在弦長(zhǎng)75%的后方，對(duì)翼型尾部進(jìn)行擺動(dòng)，如圖1所示.定義5種翼型，分別為上擺－5°、原始翼型、下擺5°、下擺10°、下擺15°.

圖1 翼型尾緣擺角及變形示意圖Fig.1 Schematic diagram of pendulum angle and deformation of airfoil

在來流方向速度U∞不變的情況下，翼型尾緣的上下擺動(dòng)改變了原始翼型的弦線，因而會(huì)對(duì)攻角αt產(chǎn)生一定變化，因此，假定尾緣向下擺動(dòng)為正，向上擺動(dòng)為負(fù)，其中，βt為尾緣相對(duì)原始翼型的擺動(dòng)角度.

在尾緣1／4處，設(shè)定擺動(dòng)角度βt變形規(guī)律為

式中，β0為原始翼型擺角，即β0＝0；βmax為尾緣最大擺角；t 為擺動(dòng)時(shí)刻.

利用二次函數(shù)來約束并形成尾緣部分的中弧線.

式中，x，y 為翼型弦長(zhǎng)c 的坐標(biāo)位置；a1，a2，a3為系數(shù).

為了得到不同擺角下的翼型，可以通過以下約束條件來確定：a.A（xA，yA）與B（xB，yB）分別為尾緣中弧線的起點(diǎn)與終點(diǎn)；b.A 點(diǎn)始終保持與中弧線相切，斜率y′A＝k；c.且kAB—＝－tan（β0＋βt）；d.保持弧長(zhǎng)不變.

根據(jù)上述條件可以得到如下方程：

式中，s0表示原始S809翼型尾緣中弧線弧長(zhǎng).

如圖2所示，假設(shè)來流為水平流動(dòng)的風(fēng)速.AB即為原弦線75%的后方，尾緣發(fā)生變形部分.OB 表示由于尾緣擺動(dòng)變形后形成的新弦線.其中，βt 即表示翼型尾緣擺動(dòng)的角度，αt表示翼型變形后的實(shí)際攻角.

圖2 翼型擺角βt與攻角αt關(guān)系簡(jiǎn)化示意圖Fig.2 Simplified schematic diagram ofβtandαt

顯然，當(dāng)βt＝0時(shí)，原弦線和新弦線重合，αt＝βt＝0.當(dāng)βt≠0時(shí)，OC＝0.75＋0.25cosβt，BC＝0.25sinβt.

由于尾緣變形而產(chǎn)生了新的弦線，攻角也隨之產(chǎn)生了變化.

最終，βt與αt之間存在著一定的函數(shù)關(guān)系，其近似表達(dá)式為

由式（8）可知，翼型尾緣擺角與攻角的變換關(guān)系如表1所示，尾緣擺角每改變5°，相應(yīng)的攻角改變約1.25°.在本文中，攻角表示來流風(fēng)速與變形產(chǎn)生的新弦線之間的夾角，向上表示正擺角，向下表示負(fù)擺角.

表1 尾緣擺角角度與攻角之間的對(duì)應(yīng)關(guān)系Tab.1 Relationship between pendulum angle and corresponding angleofattack

2 計(jì)算與結(jié)果分析

2.1 網(wǎng)格劃分與計(jì)算條件

S809原始翼型流場(chǎng)計(jì)算域如圖3所示.計(jì)算域主要由兩部分組成，分別為上游來流和下游尾跡區(qū).上游來流區(qū)是C型區(qū)域，半徑為10倍的弦長(zhǎng)；下游尾跡區(qū)是正方形區(qū)域，其邊長(zhǎng)為20倍的弦長(zhǎng).

圖3 流場(chǎng)計(jì)算域示意圖Fig.3 Schematic diagram of flow field calculation domain

網(wǎng)格質(zhì)量的優(yōu)劣對(duì)數(shù)值計(jì)算有著關(guān)鍵性的影響，對(duì)于不同的計(jì)算工況，計(jì)算域的離散方法可以分為兩類：非結(jié)構(gòu)化和結(jié)構(gòu)化網(wǎng)格.計(jì)算域整體網(wǎng)格分布如圖4（a）所示，由于在翼型前緣與尾緣附近的壓差變化較大，為捕捉翼型邊界層附近的流動(dòng)，在翼型周圍布置邊界層網(wǎng)格，如圖4（b）所示.網(wǎng)格劃分為結(jié)構(gòu)化網(wǎng)格布置，第一層網(wǎng)格尺度為0.001m，計(jì)算域網(wǎng)格總數(shù)為50 000.

圖4 計(jì)算域網(wǎng)格劃分Fig.4 Mesh of calculation domain

邊界條件的設(shè)定：a.入口邊界設(shè)定為風(fēng)速恒定的穩(wěn)態(tài)風(fēng)，湍流度選取1%；b.出口邊界采用總壓為0的壓力出口；c.壁面條件采用無滑移條件；d.湍流模型 采 用S－A（Spalart－Allmaras）模 型，雷 諾 數(shù) 取7.5×105.

2.2 可靠性驗(yàn)證

為了驗(yàn)證上述計(jì)算模型及邊界條件的可靠性，選取S809原始翼型，將其升阻力計(jì)算結(jié)果與實(shí)驗(yàn)值進(jìn)行比較.本文中采用了OSU（Ohio State University）實(shí)驗(yàn)［16］的雷諾數(shù)為7.5×105時(shí)的S809 翼型實(shí)驗(yàn)數(shù)據(jù)作為二維靜態(tài)測(cè)試數(shù)據(jù).

圖5是利用上述網(wǎng)格及邊界條件，通過Fluent模擬出的原始翼型的靜態(tài)升力系數(shù)CL與阻力系數(shù)CD在不同攻角（－6°～18°）下的變化趨勢(shì).

圖5 S809實(shí)驗(yàn)值與模擬值的對(duì)比Fig.5 Comparison of experiment and simulation result of S809

圖5（a）為Re＝7.5×105情況下的OSU 風(fēng)洞測(cè)試靜態(tài)S809翼型升力系數(shù).從圖5（a）中可以看出，最大的升力系數(shù)在約16°攻角處獲得，數(shù)值為1.02.當(dāng)攻角處于－6°～11°時(shí)，模擬的升力系數(shù)比實(shí)驗(yàn)值大4.2%；而當(dāng)攻角在11°～18°時(shí)，計(jì)算得出的CL誤差值增至9%.

從圖5（b）可以得到，實(shí)驗(yàn)值阻力系數(shù)在較大的攻角范圍內(nèi)均保持在較低的狀態(tài).當(dāng)攻角處于－6°～16°時(shí)，模擬的阻力系數(shù)與實(shí)驗(yàn)值的誤差均保持在8%以內(nèi)，當(dāng)攻角處于16°～18°時(shí)，誤差值達(dá)到10%.

由于網(wǎng)格對(duì)CFD（computational fluid dynamics）計(jì)算結(jié)果也有著至關(guān)重要的影響，故作者對(duì)網(wǎng)格數(shù)量作了驗(yàn)證，分別試用了網(wǎng)格數(shù)為34 520的非結(jié)構(gòu)性網(wǎng)格和網(wǎng)格數(shù)為38 580的結(jié)構(gòu)性網(wǎng)格對(duì)S809原始翼型進(jìn)行模擬，計(jì)算結(jié)果與實(shí)驗(yàn)值的平均誤差超過10%，即證明了50 000以下不同數(shù)量的網(wǎng)格對(duì)計(jì)算結(jié)果有著不可忽視的影響.但是，一旦網(wǎng)格數(shù)超過50 000，繼續(xù)對(duì)網(wǎng)格進(jìn)行加密，無論增加多少網(wǎng)格，對(duì)計(jì)算結(jié)果沒有太大的影響.另外，考慮計(jì)算經(jīng)濟(jì)性以及時(shí)間成本的情況下，50 000的網(wǎng)格數(shù)即滿足計(jì)算要求.

此外，在CFD計(jì)算中使用了壁面函數(shù)，因?yàn)椋诒诿嬷車膲毫μ荻群艽?，需要排除在邊界層中存在解析解的情況，并檢查了一次層網(wǎng)格質(zhì)心到壁面的無量綱距離y＋.從圖6 可以看出，y＋的范圍較小，第一層網(wǎng)格布置得比較合理.因此，采用上述網(wǎng)格數(shù)量及質(zhì)量可以得到較為準(zhǔn)確的模擬結(jié)果.

圖6 壁面附近（S809原始翼型）Fig.6 y＋near the airfoil（original S809airfoil）

綜上所述，升阻力系數(shù)的模擬結(jié)果與試驗(yàn)值曲線吻合得比較好，雖然存在一定誤差，但是，研究重點(diǎn)是不同尾緣擺角對(duì)翼型氣動(dòng)性能的影響，屬于比較性研究，稍有偏差并不影響最終結(jié)論.所以，采用的算法是可靠的.

2.3 尾緣擺角對(duì)翼型升阻力系數(shù)的影響

圖7為S809原始翼型及上擺－5°、下擺5°、下擺10°、下擺15°后翼型在不同攻角下的靜態(tài)CL和CD.

圖7 尾緣擺動(dòng)后靜態(tài)升阻力系數(shù)隨攻角的變化關(guān)系Fig.7 Variation trend of lift and drag coefficient of deformed airfoils with the change of angle of attack

由圖7（a）可知，在－6°～12°攻角范圍內(nèi)，下擺翼型的CL整體上較原始翼型高.表2 為擺角翼型相對(duì)于S809原始翼型升力系數(shù)提高的平均值.

表2 擺角翼型相對(duì)于S809原始翼型升力系數(shù)提高的平均值Tab.2 Comparison of the average value of the increasing CLof deformed airfoils with that of original S809airfoil

原始翼型的升力系數(shù)在小攻角范圍內(nèi)，會(huì)隨著攻角的變大而保持變大趨勢(shì)，其中，對(duì)于下擺5°，10°，15°這3條升力曲線，當(dāng)攻角在－5°～4°區(qū)間為線性區(qū)域（翼型表面為完全附著流動(dòng)），而攻角在4°～12°區(qū)間則為非線性區(qū)域（翼型表面為轉(zhuǎn)捩流動(dòng)或部分分離流動(dòng)）.

原始翼型在－6°～16°攻角范圍，升力系數(shù)變化范圍為－0.48～1.14，而下擺翼型當(dāng)攻角范圍為－6°～4°時(shí)，升力系數(shù)變化范圍可完全滿足原始翼型運(yùn)行于－6°～16°攻角范圍時(shí)的升力系數(shù)變化范圍，且在整個(gè)區(qū)間內(nèi)與攻角均保持線性關(guān)系（翼型表面為完全附著流動(dòng)）.

當(dāng)翼型處于5°～20°攻角范圍內(nèi)，雖然曲線呈現(xiàn)非線性特征，但是，升力系數(shù)仍隨著攻角的增大而增大.

由圖7（b）可知，對(duì)于原始翼型，當(dāng)攻角為－6°～12°時(shí)，最大與最小阻力系數(shù)變化區(qū)間為0.013 08～0.035 68；而對(duì)于尾緣擺動(dòng)后的翼型，攻角僅為－5°～4°時(shí)，最大與最小阻力系變化范圍就可以達(dá)到0.011 12～0.027 37，即在同樣升力系數(shù)區(qū)間內(nèi)，尾緣下擺后翼型的阻力系數(shù)比原始翼型的阻力系數(shù)更小.

2.4 變形翼型表面靜壓分布

如圖8所示，隨著尾緣下擺角度的增加，翼型下表面壓力也隨之增大，并且在前緣點(diǎn)和尾緣點(diǎn)位置壓力梯度特別大.從壓差的量級(jí)上來看，上擺－5°翼型上下表面壓差最?。▓D8（a）），下擺15°翼型壓差最大（圖8（e）），從而可以更清楚地證明了尾緣擺角可以提高升阻比.

圖8 變形翼型表面靜壓分布Fig.8 Static pressure distribution on the deformed airfoil surface

2.5 尾緣擺角對(duì)翼型流動(dòng)分離特性的影響

尾緣下擺對(duì)翼型升阻力系數(shù)起到了較好的改善，但是，尾緣擺動(dòng)對(duì)翼型流線存在不利的影響，如圖9所示.由圖9（a）和圖9（b）中可以看出，尾緣上擺－5°和原始翼型在8°攻角下并沒有產(chǎn)生旋渦，而下擺翼型出現(xiàn)了旋渦，對(duì)應(yīng)到升阻力系數(shù)圖（圖7）上，可以看到8°攻角處上擺翼型的升阻力系數(shù)明顯發(fā)生了突變.

圖9 翼型流線圖（8°攻角）Fig.9 Streamline of airfoil（angle of attack is 8°）

下擺5°時(shí)（圖9（c））開始出現(xiàn)較小的旋渦，并且旋渦在翼型尾部剛產(chǎn)生；下擺10°時(shí)（圖9（d）），旋渦變大，并且開始向翼型前緣處蔓延，并可以看出旋渦呈脫落狀；到下擺15°時(shí)（圖9（e）），明顯可以看到有一對(duì)旋轉(zhuǎn)方向相反的旋渦存在.從圖9（a）—9（e）的變化趨勢(shì)可以看出，隨著尾緣下擺角度的增加，翼型尾緣部分出現(xiàn)流線分離的攻角明顯變小.一旦翼型尾緣出現(xiàn)流線分離后，由于湍流模型等因素，模擬值的誤差會(huì)增大.

上擺－5°翼型相對(duì)于原始翼型，CL平均減小約17.6%，CD變化較小.另?yè)?jù)計(jì)算，原始翼型于16°攻角開始產(chǎn)生流線分離現(xiàn)象，而上擺－5°翼型在攻角14°時(shí)即開始發(fā)生流線分離.

綜上所述，原始S809 翼型高升阻比區(qū)間為－6°～16°區(qū)間，隨著S809尾緣下擺5°，10°和15°，翼型的升阻力系數(shù)較原始翼型的有較大的改善，而且最大升阻比明顯變大.然而，由于尾緣擺動(dòng)，改變了原始翼型的流動(dòng)性能，翼型隨著下擺角度的增大，翼型產(chǎn)生分離渦的攻角卻隨之提前，翼型完全附著流動(dòng)攻角區(qū)間減小到－5°～4°，正常運(yùn)行區(qū)間大幅減小，更容易形成失速.對(duì)于上擺－5°翼型，其升阻比較原始翼型的差，翼型完全附著流動(dòng)攻角區(qū)間為－3°～14°，總體而言氣動(dòng)特性較原始翼型的差.

3 結(jié) 論

對(duì)S809翼型的尾緣部分進(jìn)行適當(dāng)變形，利用CFD軟件對(duì)S809上擺－5°、下擺5°，10°及15°變形翼型作了數(shù)值模擬，分析了各個(gè)翼型的特點(diǎn)，針對(duì)翼型的升阻力系數(shù)、上下表面壓力、失速等方面，研究了攻角與擺角對(duì)翼型氣動(dòng)特性的影響，并總結(jié)了相關(guān)規(guī)律.

a.在小攻角工況下，尾緣擺角較大的翼型相對(duì)于擺角較小的翼型升力系數(shù)更高，而阻力系數(shù)基本不變，故升阻比得到了較好的提高，升力系數(shù)的最大值變大，大升阻比攻角范圍比原始翼型的更寬.

b.對(duì)于變形翼型，隨著翼型尾緣擺角的增大，其尾緣部分的流動(dòng)就越容易出現(xiàn)分離，而且，當(dāng)攻角超過4°，分離渦向翼型前緣蔓延的趨勢(shì)就越明顯，會(huì)導(dǎo)致較多的損失，最終引起阻力系數(shù)的增大.

c.隨著尾緣下擺角度變大，翼型上下表面壓差也隨之明顯變大.翼型下表面壓力增大，上表面壓力逐漸減小，并且在前緣點(diǎn)與尾緣點(diǎn)位置壓力梯度特別大.上擺－5°時(shí)翼型上下表面壓差最小，下擺15°時(shí)翼型上下表面壓差最大.

d.上擺翼型升阻力特性及失速特性均較原始翼型的差，故上擺翼型的氣動(dòng)性能不如S809原始翼型的.

［1］李春，葉舟，高偉，等.現(xiàn)代陸海風(fēng)力機(jī)計(jì)算與仿真［M］.上海：上海科學(xué)技術(shù)出版社，2012.

［2］The Danish Society of Engineers and the Federation of Engineers.DS 472Loads and safety of wind turbine construction［S］.Copenhagen：Dansk Standard，1992.

［3］Kooijman H J T，Lindenburg C，Winkelaar D，et al.DOWEC 6MW pre－design：aero－elastic modeling of the DOWEC 6 MW pre－design in PHATAS［R］.Petten：Energy Research Center of the Netherlands，2003.

［4］Sharma R N，Madawala U K.The concept of a smart wind turbine system［J］.Renewable Energy，2012，39（1）：403–410.

［5］Pasupulati S V，Wallace J，Dawson M.Variable length blades wind turbine［C］∥Power Engineering Society General Meeting，IEEE，2005，3：2097－2100.

［6］Barlas T K，Van Kuik G A M.Review of state of the art in smart rotor control research for wind turbines［J］.Progress in Aerospace Sciences，2010，46（1）：1－27.

［7］Mohamed M H，Wetzel K K.3D woven carbon／glass hybrid spar cap for wind turbine rotor blade［J］.Journal of Solar Energy Engineering，2006，128（4）：562－573.

［8］Rajadurai J S，Christopher T，Thanigaiyarasu G，et al.Finite element analysis with an improved failure criterion for composite wind turbine blades［J］.Forschung im Ingenieurwesen，2008，72（4）：193－207.

［9］張建，陳榴，戴韌，等.下游塔架對(duì)水平軸風(fēng)力機(jī)葉片氣動(dòng)負(fù)荷的作用［J］.上海理工大學(xué)學(xué)報(bào)，2010，32（1）：73－78.

［10］Andrew T L，Richard G J.Compliant blades for wind turbines［J］.IPENZ Transactions，1999，26（1）：7－12.

［11］Veers P S，Eisler G R，Laino D J，et al.The use of twist－coupled blades to enhance the performance of horizontal axis wind turbines［R］.New Mexico：Sandia National Laboratories，2001.

［12］Maheri A，Noroozi S，Toomer C A，et al.WTAB，a computer program for predicting the performance of horizontal axis wind turbines with adaptive blades［J］.Renewable Energy，2006，31（11）：1673－1685.

［13］Maheri A，Noroozi S，Vinney J.Application of combined analytical／FEA coupled aero－structure simulation in design of wind turbine adaptive blades［J］.Renewable Energy，2007，32（12）：2011－2018.

［14］Maheri A，Isikveren A T.Design of wind turbine passive smart blades［C］∥European Wind Energy Conference.2009.

［15］NREL’s S809 Airfoil graphic and coordinates［EB／OL］.NREL.（2009－10－21）［2013－4－20］.http：∥wind.nrel.gov／airfoils／Shapes／S809－Shape.html.

［16］Ramsay R R，Hoffman M J，Gregore K G.Effects of grit roughness and pitch oscillations on the S809airfoil［R］.Golden，Colorado：National Renewable Energy Laboratory，1995.