王海鵬，王華明

（南京航空航天大學(xué) 直升機(jī)旋翼動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210016）

0 引 言

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

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

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

1 計(jì)算方法

1.1 第二代低階面元法

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

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

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

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

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

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

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

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

1.2 自由尾跡計(jì)算

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

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

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

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

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

1.3 計(jì)算步驟與迭代方法

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

圖2 計(jì)算流程框圖Fig.2 Calculation flow chart

1.4 粘性修正

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

2 計(jì)算結(jié)果對(duì)比

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

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

為便于對(duì)比，本文還采用CFD方法求解雷諾時(shí)均N－S方程計(jì)算該旋翼氣動(dòng)特性。在CFD計(jì)算中采用全結(jié)構(gòu)化網(wǎng)格，網(wǎng)格數(shù)量約200萬，采用剪切應(yīng)力輸運(yùn)（SST）湍流模型。

圖3為以上兩種方法計(jì)算槳葉壓力分布與試驗(yàn)值的對(duì)比。從對(duì)比結(jié)果可以看出，以上兩種方法與試驗(yàn)結(jié)果有很好的一致性。

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

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

圖4 旋翼拉力系數(shù)對(duì)比Fig.4 The comparison chart of rotor lift coefficient

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

圖5 旋翼力矩系數(shù)對(duì)比Fig.5 The comparison chart of rotor moment coefficient

圖6為上述兩種方法計(jì)算拉力系數(shù)與試驗(yàn)值的相對(duì)誤差對(duì)比圖。從圖中可以看出，兩種方法計(jì)算精度相當(dāng)。

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

圖6 拉力系數(shù)相對(duì)誤差對(duì)比圖Fig.6 The comparison chart of lift coefficient relative error

圖7 計(jì)算耗時(shí)對(duì)比圖Fig.7 The comparison chart of calculation time cost

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

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

3 結(jié)論與展望

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

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

通過引入槳葉的揮舞運(yùn)動(dòng)，本文提出的計(jì)算方法可拓展至直升機(jī)前飛時(shí)旋翼的氣動(dòng)特性計(jì)算。本文提出的計(jì)算方法和旋翼氣動(dòng)噪聲分析方法相結(jié)合，可以預(yù)估旋翼的氣動(dòng)噪聲。

［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］徐國(guó)華.應(yīng)用自由尾跡分析的新型槳尖旋翼氣動(dòng)特性研究［D］.南京：南京航空航天大學(xué)，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ì)算方法［J］.航空學(xué)報(bào)，2004，25（6）：546－550.

［7］尹堅(jiān)平，胡章偉.旋翼表面非定常壓力脈動(dòng)計(jì)算的三維自由尾跡非定常面元法［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，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］仲唯貴.直升機(jī)旋翼尾渦與尾槳干擾噪聲研究［D］.南京：南京航空航天大學(xué)，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.