喻延福+楊志剛+賈青+Ivan+Dobrev

摘要：為在直升機(jī)旋翼氣動性能數(shù)值模擬時簡化建模過程、縮減計算時間，利用用戶自定義函數(shù)（User Defined Function，UDF）設(shè)計混合模型盤激勵模型和線激勵模型，并對簡單旋翼的懸停工況進(jìn)行模擬.與風(fēng)洞試驗結(jié)果的對比表明：所設(shè)計的混合模型在簡化旋翼模型的同時，能有效地計算旋翼的氣動特性，模擬旋翼懸停時的流場，具有正確性和可行性；盤激勵模型作為定常計算模型能夠較快地計算得到旋翼的氣動性能，缺點是不能體現(xiàn)每個槳葉對流場的單獨作用；所設(shè)計的線激勵模型在計算時由于所用的誘導(dǎo)速度為平均值，所以計算結(jié)果中旋翼效率比實際值偏高；通過與粒子圖像測速（Particle Image Velocimetry，PIV）測量結(jié)果對比發(fā)現(xiàn)，線激勵模型能較好地模擬出槳尖渦的分布.

關(guān)鍵詞：直升機(jī)； 旋翼； 氣動特性； 風(fēng)洞試驗； 用戶自定義函數(shù)； 粒子圖像測速； 數(shù)值模擬

中圖分類號： V211.52

文獻(xiàn)標(biāo)志碼： B

Abstract：To simplify the modeling process and reduce computation time during the numerical simulation of aerodynamic performance of helicopter rotor， the actuator disk model and the actuator line model of hybrid model are developed using User Defined Function（UDF） and a simple rotor in hovering is simulated. The results are compared with the results of wind tunnel test. It is indicated that， while the rotor model is simplified， these hybrid models are effective and feasible to calculate the aerodynamic characteristic of helicopter rotor in hovering and simulate the flow field； the actuator disk model is a steady solution model which is able to calculate the aerodynamic characteristic of rotor in a short time， but it cannot separately present the effect of every blade on the flow field； as the fact that an average induced velocity is used in the actuator line model， the calculated result of rotor efficiency is a little higher than the actual value； by the comparison between the hybrid model results and the Particle Image Velocimetry（PIV） measurement results， it was found that the actuator line model could well predict the distribution of blade tip vortex.

Key words：helicopter； rotor； aerodynamic characteristic； wind tunnel test； user defined function； particle image velocimetry； numerical simulation

0 引 言

對風(fēng)力機(jī)的風(fēng)輪或直升機(jī)旋翼等旋轉(zhuǎn)流體機(jī)械進(jìn)行數(shù)值計算時，如使用傳統(tǒng)的CFD方法，需將槳葉轉(zhuǎn)子幾何表面的邊界條件設(shè)為壁面，在壁面生成邊界層并在適當(dāng)區(qū)域加密，這對模型的網(wǎng)格質(zhì)量和數(shù)量有較高要求；如要研究槳尖渦和進(jìn)場尾跡，還需進(jìn)行非定常計算，不僅工作量大，且對計算資源的要求高.例如，在風(fēng)力發(fā)電農(nóng)場的氣動研究中，上游風(fēng)力機(jī)的槳葉轉(zhuǎn)子產(chǎn)生的尾流會對下游風(fēng)力機(jī)產(chǎn)生影響，為仿真研究這種現(xiàn)象，需要對多個風(fēng)力機(jī)同時進(jìn)行非定常數(shù)值計算，對計算資源有極高的要求，計算時間長.因此，有必要使用一種能加速計算同時可正確模擬槳葉與空氣間相互氣動作用、仿真尾流分布情況的仿真方法.

混合模型[1]是進(jìn)行風(fēng)力電廠CFD計算時應(yīng)用較多的一種方法.該模型簡化真實槳葉轉(zhuǎn)子模型實現(xiàn)加速計算.在該方法中，真實的槳葉轉(zhuǎn)子幾何模型被外形較簡單的特殊幾何模型代替，模型的網(wǎng)格數(shù)量較少且無須生成邊界層，能模擬真實的槳葉轉(zhuǎn)子與空氣的相互作用，對計算資源的要求較低.之所以稱之為混合模型，是因為在數(shù)值計算時將動量理論和葉素理論與數(shù)值模擬時控制方程的迭代求解相聯(lián)系：在每一步迭代結(jié)束時，通過理論方法計算氣流流經(jīng)槳葉轉(zhuǎn)子時所受的力，再將此力以源項的形式代入到控制方程的下一步迭代中，從而模擬槳葉對氣流的作用力.根據(jù)代替槳葉轉(zhuǎn)子的幾何模型不同，針對風(fēng)力機(jī)數(shù)值模擬的混合模型主要發(fā)展出3種形式：盤激勵模型 [2]、線激勵模型[3-4]和面激勵模型[5-9].在盤激勵模型中，用一個類似槳盤的圓盤代替槳葉轉(zhuǎn)子，對流經(jīng)圓盤的氣流施加單位體積力或單位面積力.因此力沿圓盤的方位角方向均勻分布而沿徑向有變化，故其缺點是無法體現(xiàn)每個葉片的單獨作用.在線激勵模型中，每根槳葉被一條沿槳盤半徑方向的線段或細(xì)桿代替，線段或細(xì)桿上所提供的力沿槳葉展向變化.與盤激勵模型相比，此模型能體現(xiàn)每個槳葉對氣流的單獨作用以及槳尾和槳尖渦的作用，改善流場特征的體現(xiàn).在面激勵模型中，每個槳葉被其中心面所代替，在數(shù)值計算時通過邊界條件在此中心面上提供壓力差.與前2個模型相比，面激勵模型模擬得到的尾流更接近真實情況，但在流速較大時，由于流場的分離作用，仿真結(jié)果不能與試驗結(jié)果很好吻合.[10]

上述研究主要集中在對風(fēng)力機(jī)的數(shù)值仿真方面，對于直升機(jī)旋翼這種通過旋翼旋轉(zhuǎn)對氣流加速、來流速度未知的槳葉轉(zhuǎn)子情況并未涉及.本文將上述風(fēng)力機(jī)風(fēng)輪模擬混合模型的方法應(yīng)用于直升機(jī)旋翼的數(shù)值模擬中，基于旋翼的氣動理論，借助用戶自定義函數(shù)（User Defined Function， UDF）設(shè)計一套針對旋翼氣動模擬的盤激勵模型和線激勵模型，并通過對某小型旋翼進(jìn)行懸停狀態(tài)的數(shù)值計算和風(fēng)洞試驗，分析計算與試驗結(jié)果，對所設(shè)計的模型予以驗證.對面激勵模型的流場分離現(xiàn)象和模擬槳葉振動時的局限性，本文不進(jìn)行討論.

1 混合模型設(shè)計

1.1 理論計算

在混合模型中，代替旋翼的簡單幾何模型通過對空氣施加單位體積力或單位面積力模擬槳葉與空氣的相互作用，UDF理論計算部分的主要目的是計算簡單幾何模型中每個微元上應(yīng)施加的力.

以盤激勵模型為例（見圖1），用圓盤代替旋翼后，根據(jù)微分原理，槳葉上寬度為dr的葉素被圓盤上徑向?qū)挾葹閐r的同心圓環(huán)微元代替.半徑為r處的葉素所施加的作用力被均勻分布在半徑同為r的圓環(huán)微元上.針對圖1中半徑為r處的葉素，其二維翼型的受力分析見圖2.

1.2 UDF編寫

利用CFD求解商業(yè)軟件FLUENT中的UDF功能實現(xiàn)混合模型中計算功能.在UDF中通過特定的函數(shù)得到每個網(wǎng)格單元的體積、流體速度矢量、三維坐標(biāo)等變量的值，進(jìn)而將其利用在自定義函數(shù)的計算中.模型中所用函數(shù)有x，y和z等3個方向的源項函數(shù)，計算旋翼的簡單幾何模型內(nèi)單位密度力，每步迭代后計算旋翼升力、扭矩功率以及旋翼效率的函數(shù).由于線激勵模型的非定常計算量較大，故在UDF進(jìn)行單節(jié)點試算成功之后，還需進(jìn)行并行化改寫，以滿足多節(jié)點并行計算的需要.加入UDF后的控制方程的迭代求解步驟見圖3.

利用葉素理論求解葉素的升阻力，需要在給定翼型和迎角時知道CL和CD.在UDF的編寫過程中需提前錄入所需翼型的升阻特性數(shù)據(jù).為與試驗保持一致，本算例選用NACA0015翼型，由于雷諾數(shù)較低，約為Re=0.1×106，故采用文獻(xiàn)[11]中NACA0015在雷諾數(shù)較低情況下的試驗數(shù)據(jù)，見圖4和5.在升阻特性曲線上截取若干數(shù)據(jù)點，UDF中通過插值函數(shù)進(jìn)行取值.線激勵模型中計算源項時，對于誘導(dǎo)速度的獲得，未用UDF獲取網(wǎng)格節(jié)點處轉(zhuǎn)軸方向的速度，原因是此處網(wǎng)格的局部氣流速度不能代表宏觀的旋翼上游誘導(dǎo)速度.本文利用動量理論計算旋翼的上游誘導(dǎo)速度的平均值，相當(dāng)于將誘導(dǎo)流場看作一種流速均勻的流場.這種流場情況只有當(dāng)槳葉有理想情況的扭轉(zhuǎn)時才能實現(xiàn)，而當(dāng)槳葉無扭轉(zhuǎn)時，這是一種會使得旋翼的效率比實際值偏高的簡化方法.[12]

2 數(shù)值計算

對具有2個無扭轉(zhuǎn)槳葉的小型旋翼進(jìn)行懸停狀態(tài)的數(shù)值計算.旋翼半徑R=290 mm，忽略槳轂部分的幾何結(jié)構(gòu)，槳葉在槳轂部分的無效區(qū)域徑向長度為41 mm.針對盤激勵模型和線激勵模型建立相應(yīng)的CFD模型，模型中網(wǎng)格區(qū)域分為2個部分，即旋翼區(qū)域和旋翼周邊圓柱形計算域.在盤激勵模型中，旋翼被厚度為0.01 mm的圓盤所代替，圓盤半徑為R.在線激勵模型中，旋翼區(qū)域為半徑360 mm，厚度40 mm的圓盤形旋轉(zhuǎn)區(qū)域，旋翼的2片槳葉被2個與槳葉等長（即L=R=290 mm）、截面半徑r=5 mm的圓桿代替，被圓盤包裹，見圖6.槳轂無效區(qū)半徑均設(shè)為41 mm；計算域半徑取3R=870 mm，上游長度取5R=1 450 mm，下游長度取15R=4 350 mm.在求解計算時借助網(wǎng)格滑移，使線激勵模型的旋轉(zhuǎn)區(qū)域網(wǎng)格以一定角速度繞z軸旋轉(zhuǎn)，其余計算域網(wǎng)格靜止.參考二維翼型壓力分布，升力的產(chǎn)生集中在翼型前半部分，故圓桿截面的位置位于槳葉截面翼型的前端，見圖7.由于盤激勵模型結(jié)構(gòu)簡單，所以圓盤內(nèi)和計算域均采用六面體結(jié)構(gòu)網(wǎng)格，計算域靠近圓盤下游的區(qū)域適當(dāng)加密，總網(wǎng)格數(shù)量約為165萬個.對于線激勵模型，在旋轉(zhuǎn)區(qū)域的2個圓桿內(nèi)部以及圓盤周圍計算域采用結(jié)構(gòu)網(wǎng)格，圓盤內(nèi)其余區(qū)域采用四面體非結(jié)構(gòu)網(wǎng)格，計算域下游槳尖渦區(qū)域網(wǎng)格適當(dāng)加密，總網(wǎng)格數(shù)量約為442萬個.

本算例馬赫數(shù)小于0.3，求解不可壓N-S方程.盤激勵模型進(jìn)行定常求解，湍流模型選擇標(biāo)準(zhǔn)k-ε模型.計算域入口邊界條件設(shè)為壓力入口，出口為壓力出口，其余壁面為無剪切應(yīng)力壁面.圓盤內(nèi)部網(wǎng)格區(qū)域設(shè)置源項UDF，圓盤表面面網(wǎng)格為內(nèi)部面網(wǎng)格.線激勵模型進(jìn)行非定常求解，湍流模型選擇k-ω SST模型.計算域出入口及壁面邊界條件與盤激勵模型相同，圓桿內(nèi)部網(wǎng)格區(qū)域設(shè)置源項UDF，圓盤內(nèi)的網(wǎng)格為動網(wǎng)格，以角速度ω繞旋翼轉(zhuǎn)軸運動，周圍計算域的網(wǎng)格靜止，利用網(wǎng)格滑移，圓盤面網(wǎng)格邊界條件設(shè)為交界面.動量、湍動能和耗散率的離散采用二階迎風(fēng)差分格式，壓力-速度耦合使用SIMPLE算法求解.

通過改變槳距角θ，利用2種混合模型分別模擬θ分別為4°，6°，8°，10°和12°時旋翼在轉(zhuǎn)速為勻速2 000 r/min時的情況.

3 風(fēng)洞試驗

3.1 試驗臺架布置

為驗證所設(shè)計的混合模型的正確性，利用試驗方法驗證數(shù)值計算結(jié)果.試驗在法國國立高等工程技術(shù)大學(xué)風(fēng)洞實驗室中的六分力天平上進(jìn)行.試驗風(fēng)洞為半開口式回流風(fēng)洞，試驗段長度為2 m，噴口大小為1.35 m×1.8 m.旋翼總成及支架安裝在試驗段的六分力測力天平上，使旋翼轉(zhuǎn)軸在試驗段中間位置，方向水平，離地高度700 mm，見圖8.旋翼模型尺寸與前文一致.通過更換槳轂上固定槳葉的鉸鏈改變θ.旋翼安裝在水平轉(zhuǎn)軸上，轉(zhuǎn)軸和電機(jī)之間通過圓柱軸扭矩傳感器T20WN連接，通過該傳感器測量扭矩和轉(zhuǎn)速.通過在控制系統(tǒng)上改變電機(jī)的輸入頻率改變旋翼轉(zhuǎn)速.由于旋翼尺寸較小且效率較低，旋翼產(chǎn)生的作用力不會在整個風(fēng)洞中產(chǎn)生循環(huán)氣流，故試驗時風(fēng)洞的風(fēng)機(jī)不開啟.開啟旋翼電機(jī)前先將天平y(tǒng)和z方向受力調(diào)零.開啟旋翼并達(dá)到穩(wěn)定速度后，微調(diào)圓形托盤使得y方向受力為0，即保證旋翼轉(zhuǎn)軸方向與z軸方向一致.試驗時，在相同θ下測量旋翼轉(zhuǎn)速不同時的升力、阻力矩和電機(jī)功率，然后更換鉸鏈改變θ進(jìn)行下一組測量.試驗測量θ分別為2°，4°，6°和10°的數(shù)據(jù)，但2°時旋翼產(chǎn)生的升力較小，天平測量精度不夠，故沒有記錄到最終結(jié)果.

3.2 測量

試驗利用粒子圖像測速（Particle Image Velocimetry，PIV）技術(shù)測量旋翼下游的流場分布.見圖9.Litron Nano公司的YAG激光器安裝在試驗段的上方，可產(chǎn)生200 mJ的激光脈沖.激光透過下方的柱透鏡可形成鉛垂方向的平板狀光路，調(diào)整角度使光路穿過旋翼轉(zhuǎn)軸.在試驗段的進(jìn)風(fēng)口導(dǎo)入由煙霧生成器生成的橄欖油微粒煙霧，微粒直徑為1～10 μm.拍攝所用的CCD相機(jī)為Dantec FlowSense 4M，分辨率2048 dpi×2048 dpi，鏡頭為Nikkor AF-S 105 mm f/2.8G ED IF.為獲得旋翼旋轉(zhuǎn)的位置，在槳葉上粘貼光學(xué)位置傳感器，連接到電腦的PIV控制系統(tǒng)，同步激光脈沖和相機(jī)快門，設(shè)定當(dāng)旋翼轉(zhuǎn)到某一角度時進(jìn)行拍攝.根據(jù)流場速度，連續(xù)拍攝2張照片的時間間隔為20 μs.拍攝時同一工況拍攝200組照片.PIV測量現(xiàn)場圖見圖10，可以明顯看到由于槳尖渦對煙霧微粒的離心作用形成的黑點.

4 結(jié)果與分析

4.1 旋翼氣動性能

2種模型在不同θ下的氣動性能的數(shù)值結(jié)果和真實旋翼的試驗結(jié)果見圖11和12.由圖11a可知：旋翼相同θ相同時線激勵模型計算得到的升力比盤激勵模型的結(jié)果大，且θ越大二者相差越大.由圖11b可知：在扭矩功率方面，兩者結(jié)果基本重合，從數(shù)值上看在θ=12°時偏差為4%，θ<12°時偏差小于1%.究其原因，是由于線激勵模型中誘導(dǎo)速度的計算采用平均值在整個槳盤平面上均勻分布的假設(shè)，從而使得旋翼效率比實際值偏高；而對于阻力矩的值，從圖5中可看出迎角為0～13°時，翼型的CD值變化很小，所以由于誘導(dǎo)速度的偏差對基于CD計算的扭矩功率影響很小.與試驗結(jié)果相比，混合模型得到的升力和功率在量級上和趨勢上均與試驗保持一致，數(shù)值上升力和扭矩功率試驗值略微偏大，與數(shù)值結(jié)果的誤差最大約為5%.分析原因，是由于實際試驗中旋翼在轉(zhuǎn)速2 000 r/min時總成會產(chǎn)生輕微共振，為避開此共振點，在試驗時實際采集的結(jié)果是在旋翼轉(zhuǎn)速約為2 200 r/min時的數(shù)據(jù).

對比分析旋翼懸停狀態(tài)3個無量綱因數(shù)的計算結(jié)果，見圖12.盤激勵模型與線激勵模型相比：對于大小與升力相關(guān)的拉力因數(shù)CT和懸停氣動效率FM，線激勵模型的值均偏大，其原因與前面所述升力偏大的原因相同；對于與扭矩功率相關(guān)的扭矩因數(shù)，二者偏差均很小.與試驗結(jié)果進(jìn)行對比時，試驗時實際轉(zhuǎn)速偏大帶來的影響在對結(jié)果的無量綱化時可被消除.因此可看出：當(dāng)θ不同時，混合模型方法得到的無量綱因數(shù)在量級上和趨勢上均與試驗結(jié)果相一致，但數(shù)值大小上試驗結(jié)果中的功率因數(shù)CP偏低較多，θ=10°時與數(shù)值結(jié)果偏差約25%.究其原因，可能是由于基于葉素理論進(jìn)行計算的混合模型沒有考慮實際旋翼旋轉(zhuǎn)時槳葉振動以及槳尖渦導(dǎo)致的功率損失.

4.2 流場分析

θ相同且轉(zhuǎn)速相同時，對2種混合模型仿真結(jié)果與風(fēng)洞試驗結(jié)果的流場進(jìn)行對比分析.本文主要展示θ=10°時旋翼勻速轉(zhuǎn)動使得流場穩(wěn)定后旋翼周圍的流場分布情況.

混合模型仿真結(jié)果在x=0平面的壓力云圖見圖13，槳盤在z=0位置，轉(zhuǎn)軸在y=0位置.

由圖13可以看出：由于旋翼的作用，槳盤兩側(cè)形成壓力差，從而對旋翼產(chǎn)生升力.盤激勵模型由一個圓盤代替旋翼，壓力沿方位角方向均勻分布，而線激勵模型由圓桿代替槳葉，壓力差集中在槳葉周圍，故圖13a中所示的壓力差遠(yuǎn)小于圖13b中圓桿兩側(cè)的壓力差.槳盤中心由于槳轂無效區(qū)的存在，兩側(cè)無壓力差.

雖然在盤激勵模型對空氣的作用力計算中可體現(xiàn)切向的分量，對流場沿旋轉(zhuǎn)切向方向也有加速，但由于盤激勵模型中旋翼對空氣的作用力沿方位角方向上均勻分布，所以不能像線激勵模型那樣體現(xiàn)每個槳葉單獨對空氣的作用.

z=0平面的速度云圖見圖14，可更清晰地看出盤激勵模型的這一缺陷.盤激勵模型的優(yōu)點是模型簡單，數(shù)值計算時進(jìn)行定常求解，故能在短時間內(nèi)算出旋翼的氣動性能參數(shù).當(dāng)需要短時間內(nèi)計算旋翼的氣動性能時，盤激勵模型有一定優(yōu)勢.

觀察槳尖渦，Q=1 200 s-2時線激勵模型等值面渦量圖見圖15.由此可看出：線激勵模型可模擬每根槳葉形成的槳尖渦，槳尖渦螺旋向下發(fā)展，螺旋半徑減小，這與理論相一致.流場在x=0平面渦量不變量Q的計算結(jié)果云圖見圖16.從圖16可看到槳尖渦的位置.試驗進(jìn)行PIV測量時，當(dāng)旋翼轉(zhuǎn)到同一槳葉處于鉛垂向上位置時進(jìn)行拍攝.用Dantec DynamicStudio 2.30對θ=10°，轉(zhuǎn)速分別為1 100和2 200 r/min時的200組拍攝照片進(jìn)行后處理，得到的平均速度場云圖見圖17.轉(zhuǎn)軸與x軸重合，槳葉位于x=300 mm的位置.圖17顯示旋翼下游的平均速度場分布和槳尖渦的位置，圖中從左向右數(shù)，每2個渦為一對，由旋翼旋轉(zhuǎn)1周時2個槳尖分別產(chǎn)生.其槳尖渦的中心位置與仿真結(jié)果一致，但由于數(shù)值模擬的數(shù)值耗散以及網(wǎng)格滑移交界面的影響，仿真結(jié)果中第2對及之后的槳尖渦顯示不明顯.

5 結(jié) 論

借助UDF建立能夠簡化建模流程、縮減計算時間、針對旋翼CFD仿真的盤激勵和線激勵混合模型，對簡單旋翼的懸停工作狀態(tài)進(jìn)行仿真計算，通過與風(fēng)洞試驗結(jié)果進(jìn)行對比，主要結(jié)論如下.

1）所設(shè)計的混合模型在簡化旋翼模型的同時能有效計算旋翼的氣動性能并仿真旋翼的懸停狀態(tài)，具有正確性和可行性，有進(jìn)一步開發(fā)和應(yīng)用價值.

2）作為定常計算模型，盤激勵模型能夠較快計算得到旋翼的氣動性能，缺點是不能像線激勵模型那樣體現(xiàn)旋翼每個槳葉對流場的單獨作用.

3）目前設(shè)計的線激勵模型在計算時因所用的誘導(dǎo)速度為平均值，故計算結(jié)果中旋翼效率比實際值偏高.這一點需在今后的研究中進(jìn)行改善.

4）通過與PIV測量結(jié)果進(jìn)行對比表明：線激勵模型能較好地模擬槳尖渦的分布情況，在研究槳尖渦對直升機(jī)工況的影響時有效.

參考文獻(xiàn)：

[1] VERMEER L J， SRENSEN J N， CRESPO A. Wind turbine wake aerodynamics[J]. Progress in Aerospace Sciences， 2003， 39（6/7）： 467-510. DOI： 10.1016/S0376-0421（03）00078-2.

[2] MIKKELSEN R. Actuator disc methods applied to wind turbines[D]. Technical University of Denmark， 2003.

[3] SRENSEN J N， SHEN W Z. Numerical modeling of wind turbine wakes[J]. Journal of Fluids Engineering， 2002， 124（2）： 393-399. DOI： 10.1115/1.1471361.

[4] DOBREV I， MASSOUH F. Etude d'un modèle hybride pour représenter l'écoulement à travers un rotor éolien[C]//Procédure de 17ème Congrès Franais de Mécanique. Troyes， 2005.

[5] MASSOUH F， DOBREV I， RAPIN M. Numerical simulation of wind turbine performance using a hybrid model[C]//Proceeding of the 44th AIAA Aerospace Sciences Meeting and Exhibit. Nevada， 2006： 1-10. DOI： 10.2514/6.2006-782.

[6] DOBREV I， MASSOUH F， RAPIN M. Actuator surface hybrid model[J]. Journal of Physics： Conference Series， 2007（75）： 12019. DOI： 10.1088/1742-6596/75/1/012019.

[7] DOBREV I， MASSOUH F， MAALOUF B. Lifting surface method for prediction of rotor vortex wake[J]. 19ème Congrès Franais de Mécanique， 2009： 24-28.

[8] SHEN W Z， ZHANG J H， SRENSEN J N. The actuator surface model： a new Navier-Stokes based model for rotor computations[J]. Journal of Solar Energy Engineering， 2009， 131（1）：11002. DOI： 10.1115/1.3027502.

[9] WATTERS C S， MASSON C. Recent advances in modeling of wind turbine wake vortical structure using a differential actuator disk theory[J]. Journal of Physics： Conference Series， 2007（75）： 12037. DOI： 10.1088/1742-6596/75/1/012037.

[10] MEMON A A. Développement d'un modèle de surface active pour améliorer la représentation des charges aérodynamiques sur une pale éolienne[D]. Arts et Métiers ParisTech， 2012.

[11] SHELDAHL R E， KLIMAS P C. Aerodynamic characteristics of seven symmetrical airfoil sections through 180-degree angle of attack for use in aerodynamic analysis of vertical axis wind turbines： SAND-80-2114[R]. Sandia National Labs， Albuquerque， USA， 1981.

[12] LEISHMAN J G. Principles of helicopter aerodynamics with CD extra[M]. Cambridge， Eng： Cambridge University Press， 2006.

（編輯 武曉英）