王海鵬，王華明

（南京航空航天大學 直升機旋翼動力學重點實驗室，江蘇 南京 210016）

0 引 言

拉普拉斯方程是分析不可壓、無粘、無旋流場的基本方程，用拉普拉斯方程計算飛行器的氣動載荷，最常用的數值方法是面元法。面元法經過半個世紀的發(fā)展，由最初的求解流場中無升力體表面壓力分布到現在的全機氣動特性計算，其精度和穩(wěn)定性已得到大幅度的提升［1］。20世紀80年代末，美國NASA推出了根據低階面元法開發(fā)的計算程序“VSAERO”［2］，該程序通過引入內部速度勢邊界條件，在減少計算量的同時較大地提高了計算精度，受到世界航空界的廣泛關注，被認為是第二代低階面元法發(fā)展的起始。

在面元法發(fā)展的同時，旋翼自由尾跡計算方法也在不斷進步，由最初的影響系數法［3］、時間步進松弛迭代法［4］、到現在的半隱式預估修松弛迭代法［5］，其計算穩(wěn)定性與精度也在逐步提高。由于研究側重點不同，國內外旋翼自由尾跡計算往往集中研究尾跡形狀的求解，對于旋翼槳葉多采用簡單的二階升力線［5－6］或不考慮槳葉厚度的升力面模型［4］，這種處理辦法不能完全模擬槳葉的三維幾何特性，也不能計算槳葉上下表面的壓力分布與速度分布。文獻［7－9］采用低階面元法耦合自由尾跡計算旋翼氣動特性，但其尾跡計算均采用啟動過程的時間步進法，穩(wěn)定性較差，收斂速度慢。隨著計算機技術的發(fā)展，雷諾時均N－S方程（RANS）也被運用于求解旋翼的氣動特性，但由于旋翼氣動環(huán)境的特殊性與計算方法本身的特點，計算仍需耗費大量的時間與資源。

本文采用第二代低階面元法耦合自由尾跡預估修正法計算旋翼在軸流狀態(tài)下的氣動載荷。在計算中考慮槳葉的三維幾何特性，采用矩形、三角形單位平面重構槳葉幾何外形；槳葉尾跡采用全展尾跡，其中包括三圈自由尾跡（近尾跡）和四圈遠尾跡；并利用槳葉翼型風洞試驗數據對氣動載荷進行粘性修正。通過和CFD方法進行對比分析，證明了本文所采用的計算方法的優(yōu)點。

1 計算方法

1.1 第二代低階面元法

面元法總控方程為拉普拉斯方程，其為橢圓型偏微分方程。對于三維拉普拉斯方程常利用格林定理，將其轉化為升力體表面的二維問題［10］。

第一代低階面元法通常只在升力體表面布置一種奇點，或在給定某種奇點分布下布置另一種未知強度奇點，通過物面不穿透邊界條件求解未知奇點強度分布。上述計算方法往往在非控制點處存在法向速度“穿透”，對計算結果有較大影響。只有通過加密升力體便面的網格數目，才能獲得較好的計算結果。

第二代低階面元法采用矩形、三角形單位平面重構槳葉幾何外形，槳葉面元上布置未知強度的面源與偶極子，槳葉后緣脫出偶極子尾渦（全展），槳葉和尾跡模型如圖1所示。通過上述源項組合，空間中任意點P處速度勢可表示為以下形式［10］。

其中s為物面坐標點矢量；R（s，p）為物面或尾跡點至計算點的距離；n為物面點處垂直于物面指向流域的法向量；Φ∞p為遠方自由流在計算點P點處的速度勢；S為槳葉表面；Sw為槳葉尾跡表面，δ（s）與μ（s）分別為布置于槳葉表面的源項與偶極子項。

圖1 面元法的槳葉模型與尾跡模型Fig.1 Blade panels and its wake model

上述方程采用混合邊界條件求解，具體即為槳葉外表面速度矢量不穿透（紐曼邊界條件）和槳葉內表面速度勢為常值（狄利克雷邊界條件）。如下式所示：

新邊界條件的引入使得槳葉表面的源強可直接顯式求出而無需方程求解，極大地減少了計算量；內部速度勢的引入減少面元間法向速度的“泄露”，降低了面網格分布對計算結果的敏感性，從而提高了計算精度。

尾跡強度使用后緣庫塔條件確定，其等于后緣上下面元偶極子強度之差。在已知尾跡形狀后，利用邊界條件并在槳葉表面、尾跡處劃分單位面元可將2式離散為線性方程組，求解此線性方程組可得未知源強與偶極子強度，進而反推出槳葉表面速度與槳葉氣動載荷。

1.2 自由尾跡計算

利用尾跡不受力條件，得到尾跡總控方程［6］：

上式為一常微分方程，其中p為尾跡點坐標矢量，r為尾跡渦線，V（p）為在空間P點處的速度矢量合，其中包括槳葉面元，尾跡偶極子面（尾跡渦線）和遠方來流的共同作用。本文采用半隱式預估修正法［6］求解尾跡節(jié)點位置。

在自由尾跡計算中，渦核模型的選取尤為重要，其直接決定尾跡計算的收斂。本文使用文獻［6，11］中推薦使用的渦核模型，渦核半徑表達式如下式所示。

式中rc為渦核半徑；rc0為初始渦核半徑，本文中選取0.1倍槳葉剖面弦長；a為 “Oseen”系數取1.225643；Rev為渦雷諾數；v為動粘性系數；t為尾跡歷程時間；a1值在0.01至0.1之間選?。?］。

在自由尾跡計算中，評判兩次迭代的尾跡均方根差值來判斷尾跡形狀是否收斂［5］。尾跡均方根差值計算如下式所示：

1.3 計算步驟與迭代方法

計算時首先劃分槳葉面元，根據動量葉素理論估算平均誘導速度，生成初始尾跡，并劃分尾跡面元；在此基礎上根據第二代低階面元法求解槳葉面元奇點強度，并求出該尾跡下槳葉氣動載荷；接下來根據槳葉與尾跡奇點強度，運用半隱式預估修正法計算旋翼新尾跡；繼而判斷尾跡形狀是否收斂，若滿足收斂標準便結束計算，未滿足則返回第二步繼續(xù)計算。詳細步驟如圖2所示。

圖2 計算流程框圖Fig.2 Calculation flow chart

1.4 粘性修正

由于基于拉普拉斯方程的面元法無法考慮粘性力，在計算旋翼旋轉力矩時，會造成較大誤差。針對上述情況，通常采用“附面層修正”方法來考慮粘性作用影響，該方法引入有粘／無粘耦合迭代，用表面源模型模擬粘性帶來的附面層位移厚度效應。該方法可較為準確地計算粘性作用，但增加了面元法的計算量與計算耗時?？紤]到旋翼工作時，槳葉翼型剖面攻角較小，可忽略粘性力對翼型法向載荷的影響，本文即利用面元法所計算的法向氣動載荷，根據相同雷諾數下槳葉剖面的二維翼型CFD計算數據或風洞試驗數據，查表反推出對應法向力系數下的軸向力系數，用于修正粘性影響。此法也可快速估算槳葉不同半徑處攻角。

2 計算結果對比

為驗證上述計算方法的正確性，本文選用廣泛引用的“Caradonna”旋翼作為計算算例，其翼型為NACA0012，展弦比為6，無扭轉，無尖削，其他參數參見文獻［12］。

在本算例計算中，對單片槳葉劃分672個面元：其中槳葉沿展向劃分10段；沿弦向上下表面劃分48段；沿槳葉厚度方向劃分4段。每片槳葉拖出3圈自由尾跡，4圈遠尾跡，單周尾跡中等角度劃分24個尾跡節(jié)點，共計1848個尾跡節(jié)點。

為便于對比，本文還采用CFD方法求解雷諾時均N－S方程計算該旋翼氣動特性。在CFD計算中采用全結構化網格，網格數量約200萬，采用剪切應力輸運（SST）湍流模型。

圖3為以上兩種方法計算槳葉壓力分布與試驗值的對比。從對比結果可以看出，以上兩種方法與試驗結果有很好的一致性。

圖4為以上兩種方法計算旋翼在1250r／min時，對應5°、8°和12°總距角時旋翼拉力系數。通過對比可以發(fā)現，采用面元法耦合自由尾跡和CFD的計算結果與試驗結果有較好的一致性。

圖3 槳葉剖面壓力分布對比Fig.3 The comparison chart of pressure distribution on blade sections

圖4 旋翼拉力系數對比Fig.4 The comparison chart of rotor lift coefficient

圖5為以上兩種方法計算旋翼在1250r／min時，對應5°、8°和12°總距角時旋翼力矩系數。其中包括面元法修正前與修正后的計算結果。由于面元法只能考慮旋翼力矩中由槳葉誘導阻力造成的旋轉力矩，無法考慮由槳葉摩擦造所成的旋轉力矩，無摩擦修正的計算結果明顯小于CFD計算結果。利用槳葉翼型風洞試驗數據修正后，槳葉力矩計算結果與CFD計算結果有較好的一致性。

圖5 旋翼力矩系數對比Fig.5 The comparison chart of rotor moment coefficient

圖6為上述兩種方法計算拉力系數與試驗值的相對誤差對比圖。從圖中可以看出，兩種方法計算精度相當。

圖7為上述兩種計算方法耗時對比圖，計算平臺處理器為“Xeon QC E7330”。通過對比，可見本文方所述面元法耗時僅為CFD計算耗時的1／10左右，大大節(jié)省了計算時間。

圖6 拉力系數相對誤差對比圖Fig.6 The comparison chart of lift coefficient relative error

圖7 計算耗時對比圖Fig.7 The comparison chart of calculation time cost

由于本文中尾跡模型為全展自由尾跡模型，其為從槳葉后緣拖出的渦面，而非一條集中渦線，直接比較尾跡形狀有一定困難。圖8給出了1250r／min時不同總距角的尾跡計算結果，通過對比可以發(fā)現隨著總距角的增加，尾跡徑向收縮更為劇烈，軸向移動速度也有所增加，這符合一般計算規(guī)律。

圖8 旋翼尾跡形狀（只顯示3圈自由尾跡）Fig.8 The figure of rotor wake geometry

3 結論與展望

第二代低階面元法耦合自由尾跡的計算方法能夠模擬旋翼槳葉的三維幾何特性，快速計算槳葉表面的壓力分布和旋翼的氣動載荷。

本文提出的的方法計算旋翼在軸流狀態(tài)下的氣動特性與CFD方法的計算精度相仿，但計算耗費的時間僅為CFD方法的十分之一。

通過引入槳葉的揮舞運動，本文提出的計算方法可拓展至直升機前飛時旋翼的氣動特性計算。本文提出的計算方法和旋翼氣動噪聲分析方法相結合，可以預估旋翼的氣動噪聲。

［1］ERICKSON L L.Panel methods an introduction［R］.NASA TP 2995，1990.

［2］MASKEW B.Program VSAERO theory document［R］.NASA CR 4023，1987.

［3］BLISS D B.A new approach to the free wake problem for hovering rotors［A］.Annual Forum Proceedings of the American Helicopter Society［C］.1985，463－477.

［4］徐國華.應用自由尾跡分析的新型槳尖旋翼氣動特性研究［D］.南京：南京航空航天大學，1990.

［5］BAGAI A.Contributions to the mathematical modeling of rotor flow－fields using a Pseudo－Implicit Free－Wake A－nalysis［D］.College Park：University of Maryland，1995.

［6］曾洪江，胡繼忠.一種新的自由渦尾跡計算方法［J］.航空學報，2004，25（6）：546－550.

［7］尹堅平，胡章偉.旋翼表面非定常壓力脈動計算的三維自由尾跡非定常面元法［J］.空氣動力學學報，1997，15（2）：185－191.

［8］AHMED S R，VIDJAJA V T.Unsteady panel method calculation of pressure distribution on BO105model rotor blades［J］.Journal of the American Helicopter Society，1998，43（1）：47－56.

［9］仲唯貴.直升機旋翼尾渦與尾槳干擾噪聲研究［D］.南京：南京航空航天大學，2006.

［10］KATZ J，PLOTKIN A.Low－speed aerodynamics，from wing theory to panel methods［M］. McGraw－Hill，1991.

［11］BHAGWAT M J，LEISHMAN J G.Generalized viscous vortex core models for application to free－vortex wake and aeroacoustic calculations［A］.Proceedings of the 58th Annual Forum of AHS International［C］.2002.

［12］CARADONNA F X，TUNG C.Experimental and analytical studies of a model helicopter rotor in hover［R］.NASA TM 81232，1981.