鄧 磊，喬志德，楊旭東，熊俊濤

(1.西北工業(yè)大學(xué)翼型、葉柵空氣動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，陜西西安 710072;2.機(jī)械與航空與太空工程系，加州大學(xué)歐文分校，CA 92697-3975)

0 引言

在風(fēng)力機(jī)葉片氣動(dòng)性能預(yù)測(cè)和外形設(shè)計(jì)中，展向使用翼型的空氣動(dòng)力學(xué)特性是重要的參數(shù)，其準(zhǔn)確性和完備性對(duì)葉片性能預(yù)測(cè)和設(shè)計(jì)具有重要的意義，因此可靠的翼型氣動(dòng)數(shù)據(jù)的準(zhǔn)備通常是設(shè)計(jì)開(kāi)始階段花費(fèi)時(shí)間最多也最困難的部分［1］，而且翼型氣動(dòng)數(shù)據(jù)的誤差甚至被認(rèn)為是風(fēng)輪載荷和性能預(yù)測(cè)中最大的誤差來(lái)源［2－4］。

近年來(lái)，為改進(jìn)厚翼型的氣動(dòng)性能和結(jié)構(gòu)效率，研究者提出了鈍后緣翼型的概念［4－7］，這類(lèi)翼型表現(xiàn)出氣動(dòng)和結(jié)構(gòu)效率方面的雙重優(yōu)點(diǎn)。在結(jié)構(gòu)方面，與相同最大相對(duì)厚度翼型相比，這類(lèi)翼型增加了截面面積和截面上的轉(zhuǎn)動(dòng)慣量;氣動(dòng)方面，鈍后緣增加了翼型的最大升力系數(shù)和升力線斜率，并減小了對(duì)前緣粗糙的敏感度［7－9］。這類(lèi)翼型可以使一部分壓力恢復(fù)發(fā)生于尾跡部分，因此減小了吸力面的逆壓梯度，從而減緩了附面層的分離趨勢(shì)，對(duì)光滑和粗糙翼面均提高了升力特性［10］。但是同時(shí)由于鈍后緣后的脫體渦的存在，使翼型阻力有非常大的增加，給此類(lèi)翼型的使用帶來(lái)風(fēng)險(xiǎn)。由于受風(fēng)洞實(shí)驗(yàn)阻塞度的限制，這類(lèi)翼型公開(kāi)發(fā)表的實(shí)驗(yàn)數(shù)據(jù)很少，且多為低雷諾數(shù)的實(shí)驗(yàn)數(shù)據(jù)［6－8，11］，而實(shí)際大型風(fēng)力機(jī)翼型的工作雷諾數(shù)通常在1×106～1×107之間;因此CFD數(shù)值模擬對(duì)這類(lèi)翼型的設(shè)計(jì)和使用變得尤為重要。近年來(lái)國(guó)外使用幾十萬(wàn)甚至將近百萬(wàn)網(wǎng)格數(shù)對(duì)這類(lèi)翼型進(jìn)行了CFD數(shù)值求解，但是一般使用固定的轉(zhuǎn)捩位置作為流場(chǎng)計(jì)算的輸入，沒(méi)有考慮轉(zhuǎn)捩位置的時(shí)均效應(yīng)［7，10－14］。國(guó)內(nèi)對(duì)這類(lèi)翼型的研究較少:沒(méi)有查到公開(kāi)發(fā)表的實(shí)驗(yàn)結(jié)果，多為基于商業(yè)軟件的數(shù)值求解［15－17］;且一般后緣厚度較小，在2%相對(duì)厚度左右，實(shí)際大型風(fēng)力機(jī)后緣厚度有時(shí)達(dá)到20% ～50%相對(duì)厚度。

本文基于RANS方程和線性穩(wěn)定性分析的en轉(zhuǎn)捩預(yù)測(cè)程序?qū)τ袑?shí)驗(yàn)結(jié)果的風(fēng)力機(jī)翼型進(jìn)行了氣動(dòng)性能計(jì)算。通過(guò)對(duì)DU97-Flat翼型的計(jì)算分析了耦合過(guò)程和網(wǎng)格數(shù)對(duì)結(jié)果的影響;計(jì)算了DU97-W-300翼型的氣動(dòng)性能并和風(fēng)洞實(shí)驗(yàn)結(jié)果進(jìn)行了比較。

1 流場(chǎng)分析

本文采用二維定常RANS方程作為流場(chǎng)分析程序。求解N-S方程空間離散采用格心格式的有限體積法，時(shí)間推進(jìn)采用五步Rungle-Kutta格式，引入二階、四階人工粘性項(xiàng)并使用當(dāng)?shù)貢r(shí)間步長(zhǎng)、隱式殘值光順、三重網(wǎng)格等加速收斂技術(shù)。湍流模型采用一方程的Spalart-Allmaras(S-A)模型。

2 轉(zhuǎn)捩位置計(jì)算和轉(zhuǎn)捩流模型

轉(zhuǎn)捩位置計(jì)算使用基于線性穩(wěn)定性分析的en方法，這也是目前使用最普遍的可行的轉(zhuǎn)捩預(yù)測(cè)方法。擾動(dòng)放大系數(shù)計(jì)算公式為:

其中，ξ為沿物面的弧長(zhǎng)，－αi為給定頻率下求得的空間放大率，R為當(dāng)前位置雷諾數(shù)。計(jì)算中由失穩(wěn)點(diǎn)向后計(jì)算一組頻率的放大曲線(本文計(jì)算15條)，找出臨界頻率及對(duì)應(yīng)的放大系數(shù)n達(dá)到臨界放大因子的位置，即為轉(zhuǎn)捩位置。

為消除或減弱“點(diǎn)轉(zhuǎn)捩”引入的數(shù)值擾動(dòng)，本文使用了轉(zhuǎn)捩過(guò)渡區(qū)模型。由于考慮了有層流分離的轉(zhuǎn)捩問(wèn)題，本文選用Walker的模型，這種模型對(duì)于有分離的轉(zhuǎn)捩過(guò)渡區(qū)模擬較好［18］:

其中:ReΔX為基于轉(zhuǎn)捩過(guò)渡區(qū)長(zhǎng)度ΔX的雷諾數(shù)，ReXtr為基于轉(zhuǎn)捩位置Xtr的雷諾數(shù)。

有適當(dāng)轉(zhuǎn)捩過(guò)渡區(qū)長(zhǎng)度后，渦粘性表達(dá)為:

其中:εtrans為過(guò)渡區(qū)的渦粘性系數(shù)，γtr為間歇因子，在0～1之間，ε為渦粘性系數(shù)。

3 網(wǎng)格數(shù)影響

圖1 DU97和DU97-Flat外形Fig.1 Airfoil shapes of DU97 and Du97-Flat

為考查網(wǎng)格數(shù)對(duì)計(jì)算結(jié)果的影響，并找到可以得到一定精度的氣動(dòng)特性并且盡量減少計(jì)算花費(fèi)，使用由DU97-W-300翼型通過(guò)在弦線50%往后對(duì)稱(chēng)增加厚度得度得到的DU97-Flat鈍后緣風(fēng)力機(jī)翼型作為算例進(jìn)行計(jì)算［19］。DU97-Flat翼型后緣厚度10%弦長(zhǎng)，改形翼型和基本翼型外形如圖1所示。計(jì)算狀態(tài):馬赫數(shù)0.165，雷諾數(shù)3.0 ×106，攻角 α =0°，上、下表面轉(zhuǎn)捩位置分別固定于0.343 和0.320。

此前的研究中，為方便網(wǎng)格生成，一般對(duì)有鈍后緣翼型使用O－網(wǎng)格，對(duì)尖后緣使用C－網(wǎng)格。本文由網(wǎng)格生成程序Gridgen生成C網(wǎng)格，來(lái)計(jì)算鈍后緣翼型。在網(wǎng)格生成時(shí)，將后緣處的直角稍作光順，以改善在直角處的網(wǎng)格質(zhì)量，這也是對(duì)此類(lèi)翼型網(wǎng)格生成所常用的方法［7，10，11，20］。由于本文使用 RANS 方程和轉(zhuǎn)捩預(yù)測(cè)程序耦合迭代的方法來(lái)求解，求解RANS方程得到附面層方程求解所需要的附面層外邊界的壓力或者速度分布，而附面層參數(shù)的精度直接影響到轉(zhuǎn)捩預(yù)測(cè)的精度，因此附面層內(nèi)的網(wǎng)格數(shù)在本文中有更重要的意義;垂直壁面方向的網(wǎng)格數(shù)應(yīng)滿足在附面層內(nèi)有60左右的網(wǎng)格，因此在生成網(wǎng)格時(shí)，第一層網(wǎng)格距壁面1.0×10－6，內(nèi)側(cè)60層網(wǎng)格等間距分布。

為便于比較網(wǎng)格數(shù)對(duì)計(jì)算結(jié)果的影響，使用一系列的C－網(wǎng)格來(lái)計(jì)算。表1為網(wǎng)格分布比較，圖2為352×96網(wǎng)格的分布。

表1 網(wǎng)格分布Table 1 Different computational grids

圖2 計(jì)算網(wǎng)格352×96Fig.2 Computational Grid of 352 ×96

圖3為計(jì)算的升力系數(shù)隨網(wǎng)格數(shù)的變化關(guān)系。可以看出隨網(wǎng)格數(shù)的增加，升力系數(shù)逐步趨于一個(gè)收斂值，但是在網(wǎng)格數(shù)更高時(shí)，升力系數(shù)又有所增加。在計(jì)算中，似乎使用更多的網(wǎng)格可以更準(zhǔn)確地計(jì)算氣動(dòng)性能，但是隨之而來(lái)的是計(jì)算花費(fèi)的急劇增加，特別是本文使用RANS方程和轉(zhuǎn)捩位置耦合迭代，計(jì)算花費(fèi)對(duì)網(wǎng)格數(shù)的依賴(lài)關(guān)系更加明顯，因此需要在計(jì)算精度和計(jì)算花費(fèi)中間折中。本文以下的計(jì)算使用448×112網(wǎng)格。

圖3 升力系數(shù)隨網(wǎng)格數(shù)變化Fig.3 Comparison of lift coefficient for different number of grid points

4 耦合求解

此前的鈍后緣氣動(dòng)性能計(jì)算中，一般使用某個(gè)轉(zhuǎn)捩位置計(jì)算程序，計(jì)算出轉(zhuǎn)捩位置坐標(biāo)后代入到RANS方程求解器中［7，10，20］。但是由于鈍后緣處的渦脫落使翼型表面的壓力分布表現(xiàn)出周期性的脈動(dòng)，而使用這種固定的轉(zhuǎn)捩位置沒(méi)有考慮這種周期脈動(dòng)的效應(yīng)對(duì)轉(zhuǎn)捩位置計(jì)算的影響，因此Standish K.J.等建議使用有時(shí)均性的轉(zhuǎn)捩位置模型。本文耦合二維RANS方程求解和基于線性穩(wěn)定性分析的en轉(zhuǎn)捩位置計(jì)算方法，計(jì)算鈍后緣翼型的氣動(dòng)性能。此前對(duì)定常流場(chǎng)的耦合過(guò)程中，通常是將流場(chǎng)計(jì)算至某個(gè)流動(dòng)參數(shù)相對(duì)迭代次數(shù)的導(dǎo)數(shù)為常數(shù)時(shí)，而后調(diào)用轉(zhuǎn)捩位置計(jì)算［21－22］。本文中由于鈍后緣處的流場(chǎng)的非定常，不能得到收斂的流場(chǎng)解，因此為在轉(zhuǎn)捩位置計(jì)算中考慮周期脈動(dòng)的時(shí)均性，使用下面的耦合方法:

1)首先給一個(gè)初始的轉(zhuǎn)捩位置代入RANS方程計(jì)算中。建議初始轉(zhuǎn)捩位置在距前緣位置的遠(yuǎn)下游，比如80%弦長(zhǎng)或者后緣處;

2)使用這個(gè)固定的轉(zhuǎn)捩位置迭代計(jì)算RANS方程，至最大迭代次數(shù)Cyclemax，此時(shí)流場(chǎng)已經(jīng)進(jìn)入周期脈動(dòng);

3)繼續(xù)計(jì)算流場(chǎng)Δkcyc次，并求解Δkcyc內(nèi)的平均壓力分布以及平均附面層外邊界壓力和速度分布;

4)將計(jì)算的附面層外邊界速度和壓力分布作為輸入條件，調(diào)用附面層方程的求解程序，計(jì)算得到附面層參數(shù)和速度型;

5)使用計(jì)算得到的附面層參數(shù)和速度型作為輸入條件，進(jìn)行附面層流動(dòng)的穩(wěn)定性和轉(zhuǎn)捩位置求解。在從失穩(wěn)點(diǎn)向下游進(jìn)行轉(zhuǎn)捩位置計(jì)算的同時(shí)，判斷是否出現(xiàn)了由于逆壓梯度過(guò)大產(chǎn)生的層流分離，如果發(fā)生了層流分離則不再繼續(xù)向下游進(jìn)行轉(zhuǎn)捩位置計(jì)算，而以層流分離點(diǎn)作為轉(zhuǎn)捩位置。但是，此時(shí)計(jì)算的(j表示轉(zhuǎn)捩位置迭代的次數(shù))被認(rèn)為是欠松弛的，即:實(shí)際用于RANS方程使用的轉(zhuǎn)捩位置坐標(biāo)位于的下游，用下面的方法來(lái)確定:

其中，frelax為 0.7。

這樣，就完成一次轉(zhuǎn)捩位置耦合迭代過(guò)程。

本文求解中，初始轉(zhuǎn)捩位置上、下表面均設(shè)置為0.8，流場(chǎng)迭代1500次后，計(jì)算200次流場(chǎng)的平均壓力分布;以后每200步調(diào)用一次轉(zhuǎn)捩預(yù)測(cè)程序。圖4為轉(zhuǎn)捩位置隨轉(zhuǎn)捩位置迭代次數(shù)的關(guān)系，可以看出，轉(zhuǎn)捩位置在10次迭代之內(nèi)，就可以達(dá)到收斂。

圖4 轉(zhuǎn)捩位置收斂過(guò)程Fig.4 Convergence history of the transition locations

圖5 氣動(dòng)性能收斂過(guò)程Fig.5 Convergence history of lift and drag coefficient

圖5為升力系數(shù)和阻力系數(shù)隨迭代次數(shù)的關(guān)系?？梢钥闯鲭S迭代的進(jìn)行，升力系數(shù)和阻力系數(shù)均可以收斂到一個(gè)值，同時(shí)也可以看出，氣動(dòng)特性和轉(zhuǎn)捩位置的依賴(lài)關(guān)系，也驗(yàn)證了進(jìn)行轉(zhuǎn)捩位置迭代的必要性。

5 氣動(dòng)性能計(jì)算

大雷諾數(shù)的鈍后緣翼型風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)很少，1999年Timmer等發(fā)表了一個(gè)DU97-W-300翼型雷諾數(shù)為3.0 ×106的風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)［23］。DU97-W-300 翼型最大厚度30%，后緣厚度1.74%弦長(zhǎng)，相比于其他鈍后緣翼型相對(duì)較尖。本文計(jì)算網(wǎng)格448×112，計(jì)算狀態(tài)馬赫數(shù)為0.165，雷諾數(shù) 3.0 ×106。

文獻(xiàn)［12］使用RANS方程求解器SACCARA程序，進(jìn)行了DU97-W-300翼型的計(jì)算，但是使用了將近40萬(wàn)的網(wǎng)格，而轉(zhuǎn)捩位置是使用Xfoil程序求解出轉(zhuǎn)捩位置，而后作為固定轉(zhuǎn)捩代入流場(chǎng)求解器中。本文進(jìn)行了耦合求解并和［12］文計(jì)算結(jié)果作了比較。圖6為本文計(jì)算結(jié)果和Xfoil計(jì)算的轉(zhuǎn)捩位置的比較，可以看出在攻角較小時(shí)下表面轉(zhuǎn)捩位置吻合較好，但是隨攻角增加差別逐漸增大;本文計(jì)算的上表面的轉(zhuǎn)捩位置小于Xfoil計(jì)算結(jié)果，特別是本文計(jì)算中，在12°就發(fā)生了前緣分離，而Xfoil計(jì)算顯示還有10%弦長(zhǎng)左右的層流。圖7為本文CFD計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果以及文獻(xiàn)［12］計(jì)算結(jié)果的比較?？梢钥闯觯ο禂?shù)在線性段和實(shí)驗(yàn)結(jié)果吻合很好，且比參考文獻(xiàn)結(jié)果吻合好;但是和文獻(xiàn)結(jié)果一樣，失速迎角沒(méi)有捕捉到。力矩系數(shù)可以看出，在13°之前本文計(jì)算結(jié)果和實(shí)驗(yàn)值吻合很好。為考察阻力計(jì)算的精度，圖8為計(jì)算結(jié)果、實(shí)驗(yàn)結(jié)果和文獻(xiàn)結(jié)果的極曲線的比較?？梢钥闯觯诘妥璺秶鷥?nèi)阻力系數(shù)和實(shí)驗(yàn)值吻合較好。

6 結(jié)論

本文耦合RANS方程求解和轉(zhuǎn)捩位置預(yù)測(cè)程序，進(jìn)行了大厚度鈍后緣風(fēng)力機(jī)翼型氣動(dòng)性能計(jì)算研究?？梢缘玫饺缦陆Y(jié)論:

(1)通過(guò)分析網(wǎng)格數(shù)對(duì)計(jì)算結(jié)果的影響，可以看出隨網(wǎng)格數(shù)的增加，氣動(dòng)性能逐步趨于一個(gè)近似定常的值。在實(shí)際計(jì)算中，需要在計(jì)算花費(fèi)和計(jì)算精度中間找到一個(gè)折中的網(wǎng)格數(shù)。

圖6 轉(zhuǎn)捩位置比較Fig.6 Comparison of transition locations

圖7 升力系數(shù)和力矩系數(shù)比較Fig.7 Comparison of lift and moment coefficient

圖8 極曲線比較Fig.8 Comparison of drag polar

(2)本文發(fā)展的耦合方法和考慮了時(shí)均效應(yīng)的轉(zhuǎn)捩模型可以用于鈍后緣翼型的氣動(dòng)性能計(jì)算中，在10次轉(zhuǎn)捩位置迭代之內(nèi)可以得到收斂的轉(zhuǎn)捩位置并得到翼型的氣動(dòng)性能。

(3)通過(guò)本文計(jì)算的氣動(dòng)性能和實(shí)驗(yàn)結(jié)果以及文獻(xiàn)結(jié)果的比較，表明本文的耦合求解方法可以通過(guò)遠(yuǎn)少于參考文獻(xiàn)的網(wǎng)格數(shù)得到和實(shí)驗(yàn)結(jié)果吻合更好的計(jì)算結(jié)果，但是失速迎角和最大升力系數(shù)的捕捉不好，還需要進(jìn)一步的改進(jìn)。阻力和力矩的計(jì)算結(jié)果也和實(shí)驗(yàn)值吻合較好，為此類(lèi)翼型的設(shè)計(jì)和使用提供了數(shù)值模擬基礎(chǔ)。

［1］MARSHALL L.Usage advice［EB/OL］.http://wind.nrel.gov/designcodes/advice.html，Last modified 05-July-2005;accessed 28-Janu.-2010.

［2］PATRICK J.AeroDyn theory manual［R］.NREL/EL-500-36881，2005.

［3］SIMMS D，SCHRECK S，HAND M，F(xiàn)INGERSH L J.NREL unsteady aerodynamics experiment in the NASA-Ames wind tunnel:A comparison of predictions to measurements［R］.NREL/TP-500-29494.Golden，CO:National Renewable Energy Laboratory，2001.

［4］TANGLER J L.The nebulous art of using wind-tunnel airfoil data for predicting rotor performance［R］.NREL/CP-500-31243.Golden，CO:National Renewable Energy Laboratory，2002.

［5］TPI Composites.Innovative design approaches for large wind turbine blades［R］.SAND2003-0723，2003.

［6］TPI Composites.Innovative design approaches for large wind turbine blades-final report［R］，SAND2004-0074，2004.

［7］STANDISH K J，van DAM C P.Aerodynamic analysis of blunt trailing edge airfoils［J］.Journal of Solar Energy Engineering，2003，125(4):479-487.

［8］JACKSON K，ZUTECK M，van DAM C P，STANDISH K J，BERRY D.Innovative design approaches for large wind turbine blades［J］.Wind Energy，2005，8(2):141-171.

［9］Van ROOIJ R P J O M，TIMMER W A.Roughness sensitivity considerations for thick rotor blade airfoils［J］.Journal of Solar Energy Engineering，2003，125(4):468-478.

［10］BAKER J，MAYDA E，van DAM C P.Computational and experimental analysis of thick fatback wind turbine airfoils［R］.AIAA Paper 2006-193，2006.

［11］CHAO D D，van DAM C P.RaNS analysis of an inboard flatback modification on the NREL phase VI rotor［R］.AIAA Paper 2006-195，2006.

［12］WINNEMOLLER T，van DAM C P.Design and numerical optimization of thick airfoils［R］.AIAA 2006-238，2006.

［13］TPI Composites.Computational design and analysis of flatback airfoil wind tunnel experiment［R］.SAND 2008-1782，2008.

［14］STONE C，BARONE M，LYNCH C E，et al.A computational study of the aerodynamics and aeroacoustics of a flatback airfoil using hybrid RANS-LES［R］，AIAA 2009-273，2009.

［15］劉雄，陳嚴(yán)，葉枝全.增加風(fēng)力機(jī)葉片翼型后緣厚度對(duì)氣動(dòng)性能的影響［J］.太陽(yáng)能學(xué)報(bào)，2007，27(5):489-495.

［16］李秋悅，申振華.翼型進(jìn)行鈍后緣修改后氣動(dòng)性能的數(shù)值研究［J］.沈陽(yáng)航空工業(yè)學(xué)院學(xué)報(bào)，2007，24(1):1-5.

［17］張磊，楊科，趙曉路，等.不同后緣改型方式對(duì)風(fēng)力機(jī)鈍后緣翼型氣動(dòng)性能的影響［J］.工程熱物理學(xué)報(bào)，2009，30(5):773-776.

［18］HANS W S，HASSE W.Navier-Stokes airfoil computations with transition prediction including transitional flow regions［J］.AIAA Journal，2000:38(11):2059-2066.

［19］MATTHEW F B，DALE B.Aerodynamic and aeroacoustic properties of a flatback airfoil:an update［R］.AIAA Paper 2009-271，2009.

［20］STANDISH K J，van DAM C P.Analysis of blunt trailing edge airfoils［R］.AIAA 2003-353，2003.

［21］KRUMBEIN A.Navier-Stokes airfoil computation with automatic transition prediction using the DLR TAU Code-A sensitivity study［M］.New Results in Numerical and Experimental Fluid Mechanics V.，Springer Berlin/Heidelberg，Vol.92，2006.

［22］Der LEE J，JAMESON A.Natural-laminar-flow airfoil and wing design by adjoint method and automatic transition prediction［R］.AIAA 2009-897，2009.

［23］TIMMER W A，van ROOIJ R P J O M.Design and wind tunnel test of airfoil DU 97-W-300［R］.IW-98003R，Institute for Wind Energy，the Netherlands:Delft University of Technology，1999.