基于用戶自定義函數的直升機旋翼CFD模擬加速方法的設計和驗證

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

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

關鍵詞：直升機； 旋翼； 氣動特性； 風洞試驗； 用戶自定義函數； 粒子圖像測速； 數值模擬

中圖分類號： V211.52

文獻標志碼： 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 引 言

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

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

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

1 混合模型設計

1.1 理論計算

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

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

1.2 UDF編寫

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

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

2 數值計算

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

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

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

3 風洞試驗

3.1 試驗臺架布置

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

3.2 測量

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

4 結果與分析

4.1 旋翼氣動性能

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

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

4.2 流場分析

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

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

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

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

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

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

5 結 論

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

1）所設計的混合模型在簡化旋翼模型的同時能有效計算旋翼的氣動性能并仿真旋翼的懸停狀態，具有正確性和可行性，有進一步開發和應用價值.

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

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

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

參考文獻：

[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.

（編輯 武曉英）