王正之,杜孟華,朱春玲,劉森云

(1.南京工程學(xué)院能源與動(dòng)力工程學(xué)院，南京，211167；2.南京航空航天大學(xué)航空學(xué)院，南京，210016；3.中國(guó)空氣動(dòng)力研究與發(fā)展中心結(jié)冰與防除冰重點(diǎn)實(shí)驗(yàn)室，四川綿陽(yáng)，621000)

直升機(jī)相較于固定翼飛機(jī)而言飛行高度低，前飛的飛行速度慢，非常容易遭遇結(jié)冰氣象環(huán)境。直升機(jī)旋翼槳葉結(jié)冰后，改變了槳葉翼型剖面的形狀，旋翼的空氣動(dòng)力學(xué)特性被破壞，使得升力降低，阻力增加。同時(shí)，結(jié)冰后，旋翼在高速旋轉(zhuǎn)，每個(gè)截面提供的升力被改變，使升力分布不再對(duì)稱(chēng)，加劇了直升機(jī)的振動(dòng)，危及直升機(jī)的安全。報(bào)道表明我國(guó)現(xiàn)役的直升機(jī)中，諸如米-8、超黃蜂等在冬季飛行時(shí)，極易遭遇結(jié)冰氣象條件，導(dǎo)致結(jié)冰事故[1-2]。因此結(jié)冰嚴(yán)重威脅到了直升機(jī)的飛行安全，需要引起高度重視。

直升機(jī)旋翼結(jié)冰過(guò)程相對(duì)于固定翼結(jié)冰而言非常復(fù)雜。隨著計(jì)算流體力學(xué)技術(shù)的逐漸發(fā)展，相關(guān)理論的完善和成熟，數(shù)值模擬方法模擬旋翼槳葉結(jié)冰問(wèn)題并分析其對(duì)氣動(dòng)特性的影響成為了可能。在結(jié)冰對(duì)直升機(jī)氣動(dòng)特性影響數(shù)值模擬方面，美國(guó)NASA于2012年開(kāi)發(fā)的LEWICE軟件增加了直升機(jī)旋翼結(jié)冰模擬模塊[3]，首先模擬三維旋翼流場(chǎng)，然后截取其中的二維翼型流場(chǎng)結(jié)果，在此基礎(chǔ)上采用固定翼翼型結(jié)冰數(shù)值計(jì)算方法預(yù)測(cè)冰形，最后將多個(gè)截面的計(jì)算結(jié)果組合并拓展至三維，利用三維結(jié)冰外形求解結(jié)冰后的氣動(dòng)特性。目前該方法已經(jīng)應(yīng)用在UH-1H直升機(jī)懸停的試驗(yàn)程序中。在國(guó)內(nèi)，胡立芃[4]計(jì)算了旋翼槳葉不同位置翼型的結(jié)冰后氣動(dòng)特性變化，基于動(dòng)量-葉素理論計(jì)算了旋翼氣動(dòng)特性，分析了結(jié)冰對(duì)直升機(jī)飛行性能的影響。胡國(guó)才等[5]分析了槳葉結(jié)冰對(duì)尾槳?dú)鈩?dòng)特性的影響，根據(jù)某直升機(jī)尾槳槳葉結(jié)冰實(shí)驗(yàn)結(jié)果，分析了槳葉翼型結(jié)冰前后氣動(dòng)特性，利用動(dòng)量葉素理論計(jì)算了尾槳結(jié)冰前后的氣動(dòng)特性，結(jié)果表明結(jié)冰使得升力下降阻力上升，拉力系數(shù)下降而扭矩上升，氣動(dòng)性能惡化。李國(guó)知、曹義華等[6-7]以旋翼翼型冰風(fēng)洞試驗(yàn)結(jié)果為基礎(chǔ)，考慮了槳葉結(jié)冰脫落，研究了結(jié)冰對(duì)翼型氣動(dòng)特性的影響。 基于動(dòng)量-葉素理論和渦流理論，以UH-60A為模型建立了旋翼結(jié)冰后的直升機(jī)飛行動(dòng)力學(xué)模型，研究了直升機(jī)結(jié)冰后的平衡特性和穩(wěn)定性的影響，分析了結(jié)冰時(shí)間、溫度、LWC、MVD對(duì)旋翼結(jié)冰后的直升機(jī)操縱特性、姿態(tài)敏捷性、軸間耦合特性以及垂直軸操縱功效的影響。楊虎[8]以旋翼翼型冰風(fēng)洞的試驗(yàn)數(shù)據(jù)作為依據(jù)，研究了結(jié)冰對(duì)飛行力學(xué)的影響，對(duì)直升機(jī)在旋翼結(jié)冰后的平衡性、穩(wěn)定性進(jìn)行了綜合性的分析。劉國(guó)強(qiáng)等[9]分析了旋翼結(jié)冰對(duì)黑鷹直升機(jī)飛行性能影響，基于西科斯基公司翼型冰風(fēng)洞結(jié)冰試驗(yàn)推導(dǎo)的工程公式，計(jì)算了翼型結(jié)冰后的氣動(dòng)特性；在此基礎(chǔ)上建立了黑鷹直升機(jī)懸停和前飛計(jì)算模型。同時(shí)分析了槳葉結(jié)冰隨飛行速度、溫度和LWC變化規(guī)律，沿槳葉展向結(jié)冰厚度的變化情況，以及槳葉不同位置處空氣動(dòng)力對(duì)結(jié)冰的抑制作用[10]?？梢钥闯觯壳爸鄙龣C(jī)旋翼槳葉結(jié)冰后氣動(dòng)特性數(shù)值計(jì)算主要是以二維翼型計(jì)算為主，由于目前數(shù)值方法較難直接獲得槳葉三維結(jié)冰外形，較少針對(duì)三維槳葉直接開(kāi)展結(jié)冰后氣動(dòng)特性分析。

本文在前期三維旋翼結(jié)冰數(shù)值模擬研究的基礎(chǔ)之上，采用CFD方法建立了三維旋翼槳葉結(jié)冰數(shù)值模擬方法，針對(duì)結(jié)冰后的三維幾何計(jì)算了結(jié)冰前后槳葉的氣動(dòng)特性，分析了結(jié)冰對(duì)直升機(jī)旋翼槳葉升阻特性以及懸停效率等的影響。

1 計(jì)算模型

對(duì)于旋翼槳葉流場(chǎng)的計(jì)算，由于槳葉相對(duì)于地面坐標(biāo)系是旋轉(zhuǎn)運(yùn)動(dòng)的，如果直接以地面坐標(biāo)系為參考系進(jìn)行計(jì)算，那么槳葉的運(yùn)動(dòng)是一個(gè)非穩(wěn)態(tài)的運(yùn)動(dòng)，這不僅涉及到動(dòng)網(wǎng)格的問(wèn)題，而且計(jì)算量和難度都很大。

為了解決旋翼槳葉旋轉(zhuǎn)問(wèn)題，本文采用多參考系模型(MRF)對(duì)旋翼流場(chǎng)進(jìn)行求解。MRF模型將計(jì)算區(qū)域幾何劃分為多個(gè)子區(qū)域，每個(gè)子區(qū)域有自己的運(yùn)動(dòng)方式，求解區(qū)域之間信息通過(guò)交界面進(jìn)行傳遞。將整個(gè)區(qū)域劃分為2個(gè)部分，將包圍槳葉的內(nèi)場(chǎng)區(qū)域定義為旋轉(zhuǎn)坐標(biāo)系區(qū)域，遠(yuǎn)離葉片的遠(yuǎn)場(chǎng)區(qū)域?yàn)榈孛孀鴺?biāo)系。將坐標(biāo)系定義與槳葉一同旋轉(zhuǎn)則槳葉周?chē)牧鲌?chǎng)將是一個(gè)定常狀態(tài)。同時(shí)考慮到遠(yuǎn)場(chǎng)區(qū)域的流場(chǎng)是靜止?fàn)顟B(tài)，所以將遠(yuǎn)場(chǎng)區(qū)域定義為地面坐標(biāo)系。整個(gè)計(jì)算域示意如圖 1所示。

圖1 旋翼槳葉流場(chǎng)計(jì)算域

2 結(jié)冰數(shù)值計(jì)算方法

三維旋翼結(jié)冰數(shù)值模擬流程包括網(wǎng)格劃分、空氣流場(chǎng)計(jì)算、水滴撞擊特性計(jì)算、結(jié)冰生成計(jì)算和幾何模型重構(gòu)?；诘玫降娜S幾何外形就可以計(jì)算結(jié)冰前后旋翼槳葉的氣動(dòng)特性。

2.1 空氣流場(chǎng)計(jì)算

采用旋轉(zhuǎn)坐標(biāo)系下以絕對(duì)速度表示的N-S方程計(jì)算旋翼空氣流場(chǎng)，方程可以表示為以下形式：

(1)

(2)

(3)

2.2 水滴撞擊特性計(jì)算

基于空氣流場(chǎng)計(jì)算的結(jié)果，本文采用歐拉法計(jì)算水滴在旋轉(zhuǎn)過(guò)程中的運(yùn)動(dòng)軌跡，對(duì)水滴列出質(zhì)量，動(dòng)量方程，可以得到如下方程：

(4)

Kαρ(ua-uw)+αρF

(5)

式(4)～(5)中：ρw為水滴密度；uw為水滴速度；F為所受外力，這里表示重力加速度；α為水滴容積分?jǐn)?shù)；K為慣性因子。

水收集系數(shù)β是表征水滴撞擊特性的重要指標(biāo)，其定義為：微元上實(shí)際收集到的水量與可能收集到的最大水量之比，基于歐拉方法的水滴收集系數(shù)定義為：

(6)

式中：αn表示當(dāng)?shù)氐乃稳莘e分?jǐn)?shù)；α表示來(lái)流的水滴容積分?jǐn)?shù)；V和n分別表示當(dāng)?shù)氐乃嗡俣仁噶亢捅砻娣ㄏ嗍噶浚籚表示水滴來(lái)流速度，對(duì)于旋翼結(jié)冰計(jì)算來(lái)說(shuō)，水滴是相對(duì)旋翼槳葉運(yùn)動(dòng)進(jìn)而撞擊在槳葉表面，其來(lái)流速度應(yīng)由轉(zhuǎn)速來(lái)表征，因此可以采用槳尖速度表示水滴來(lái)流速度。

2.3 結(jié)冰生成計(jì)算

采用Messinger控制容積方法，將表面劃分為多個(gè)控制體，對(duì)表面每個(gè)控制體列出質(zhì)量和能量守恒方程：

(7)

(8)

2.4 幾何模型重構(gòu)

在數(shù)值模擬中網(wǎng)格劃分過(guò)程會(huì)將原本的幾何模型離散成眾多的網(wǎng)格點(diǎn)，在此基礎(chǔ)之上開(kāi)展計(jì)算。在結(jié)冰數(shù)值模擬中，如何將結(jié)冰后的外形重構(gòu)為完整幾何成為計(jì)算結(jié)冰后氣動(dòng)特性的重要一環(huán)。

根據(jù)結(jié)冰后三維槳葉表面坐標(biāo)及其拓?fù)潢P(guān)系，通過(guò)網(wǎng)格相鄰單元排布將網(wǎng)格節(jié)點(diǎn)按照弦向方向和展向方向排列，采用雙三次非均勻B樣條曲面重構(gòu)算法，重新構(gòu)建新的三維槳葉幾何外形[11]。隨后利用得到的幾何模型計(jì)算結(jié)冰后旋翼槳葉的氣動(dòng)特性。

本文研究重點(diǎn)在于分析結(jié)冰前后旋翼槳葉的氣動(dòng)特性變化情況，對(duì)旋翼結(jié)冰數(shù)值模擬方法的詳細(xì)介紹和模擬方法正確性驗(yàn)證可以參考作者前期研究成果[12-13]。

3 計(jì)算域和網(wǎng)格劃分

計(jì)算采用的旋翼槳葉為Caradonna-Tung(C-T)槳葉，槳葉的展長(zhǎng)為1.143 m，弦長(zhǎng)為0.190 5 m，展弦比為6，槳葉剖面形狀為NACA0012翼型。流場(chǎng)的邊界需要對(duì)槳葉的流場(chǎng)沒(méi)有影響，所以計(jì)算域的出入口要與槳葉保持一定的距離。將流場(chǎng)計(jì)算域入口定義在距槳葉前緣20倍旋翼槳葉展長(zhǎng)處，而流場(chǎng)計(jì)算域出口邊界定義在后緣15倍展長(zhǎng)處，開(kāi)口邊界距翼型為20倍展長(zhǎng)。在模擬過(guò)程中，定義入口邊界的氣流速度，出口邊界為壓力出口，其余邊界可使氣流自由的進(jìn)出?？紤]到C-T旋翼有2片槳葉，為了減少網(wǎng)格量以及計(jì)算量，采用了周期性邊界條件。

采用結(jié)構(gòu)網(wǎng)格對(duì)槳葉進(jìn)行劃分，生成的計(jì)算網(wǎng)格如圖 2所示。為了更好地計(jì)算槳葉的氣動(dòng)參數(shù)，對(duì)近壁面網(wǎng)格進(jìn)行了加密。為了驗(yàn)證網(wǎng)格無(wú)關(guān)性，以未結(jié)冰旋翼槳葉為例，采用不同數(shù)量的網(wǎng)格開(kāi)展流場(chǎng)計(jì)算?？偣策x取了3 240 000，5 520 000和6 240 000這3種數(shù)量的網(wǎng)格，從表1可以看出，隨著網(wǎng)格量的逐漸增多，總扭矩逐漸降低，這說(shuō)明網(wǎng)格量對(duì)計(jì)算結(jié)果有很大影響。5 520 000和6 240 000網(wǎng)格計(jì)算出的旋翼槳葉的總扭矩分別為47.8 N·m和47.6 N·m，總體結(jié)果相差不大。為了減少計(jì)算量，同時(shí)保證計(jì)算結(jié)果的準(zhǔn)確性，選擇5 520 000網(wǎng)格量來(lái)繪制網(wǎng)格。保持結(jié)冰槳葉和未結(jié)冰槳葉網(wǎng)格節(jié)點(diǎn)分布保持一致，因此在進(jìn)行結(jié)冰前后槳葉的氣動(dòng)性能對(duì)比時(shí)可以消除網(wǎng)格量的影響。

表1 網(wǎng)格無(wú)關(guān)性驗(yàn)證

圖2 旋翼槳葉網(wǎng)格劃分

4 計(jì)算結(jié)果及分析

本節(jié)利用數(shù)值模擬方法得到結(jié)冰后三維旋翼幾何外形，計(jì)算了結(jié)冰前后旋翼空氣流場(chǎng)，對(duì)比分析了結(jié)冰前后旋翼氣動(dòng)性能。旋翼轉(zhuǎn)速Ω=1 500 r/m，槳尖馬赫數(shù)Matip=0.520，槳距角為8°。

整體氣動(dòng)性能對(duì)比如下。圖3給出了懸停狀態(tài)下結(jié)冰前后旋翼槳葉上表面壓力分布云圖。從圖中可以看出，結(jié)冰后槳葉前緣出現(xiàn)了較為明顯的負(fù)壓區(qū)域，這是由于旋翼結(jié)冰破壞了槳葉原有的氣動(dòng)外形，結(jié)冰外形遮擋改變了氣流流動(dòng)情況，導(dǎo)致冰形后方出現(xiàn)了較為明顯的負(fù)壓區(qū)域，負(fù)壓區(qū)域的出現(xiàn)會(huì)對(duì)旋翼升阻力產(chǎn)生較大影響。

圖3 結(jié)冰前后槳葉表面壓力分布的影響對(duì)比

表2中給出了結(jié)冰前后旋翼的氣動(dòng)特征參數(shù)。從表中可以明顯看出結(jié)冰后旋翼的升力減小了8.8%，扭矩增加了11.1%，升力減小明顯，扭矩急劇增大。衡量旋翼系統(tǒng)氣動(dòng)效率的重要參數(shù)懸停效率從0.670 8降低到0.525 7，降幅達(dá)21.6%，說(shuō)明結(jié)冰對(duì)旋翼氣動(dòng)性能產(chǎn)生較大的影響，會(huì)導(dǎo)致直升機(jī)整體性能大幅下降。本文在計(jì)算過(guò)程中假設(shè)結(jié)冰前后旋翼實(shí)度不變，如果考慮結(jié)冰后旋翼實(shí)度增加，那么旋翼氣動(dòng)特性下降將更為劇烈。

表2 結(jié)冰對(duì)旋翼懸停性能的影響

圖4為旋翼槳葉各截面拉力系數(shù)沿展向長(zhǎng)度的分布情況，可以看出，相比于干凈旋翼，結(jié)冰后槳葉的拉力系數(shù)出現(xiàn)了下降，特別是在槳葉尖部下降幅度較大，說(shuō)明結(jié)冰對(duì)各個(gè)截面的拉力都產(chǎn)生了影響。

圖4 不同位置槳葉拉力系數(shù)沿展向分布的對(duì)比

圖5給出了不同截面處槳葉壓力系數(shù)分布曲線(xiàn)，截面位置分別為r/R=0.70,0.80,0.91,0.96。從圖中可以看出，結(jié)冰前后壓力系數(shù)分布發(fā)生較明顯的變化，結(jié)冰后壓力系數(shù)的最大值雖然增加，但是在壓力系數(shù)峰值后方壓力迅速減小，導(dǎo)致此處出現(xiàn)了明顯的壓差阻力，使得旋翼阻力大幅增加，而整個(gè)壓力系數(shù)曲線(xiàn)包圍的面積也有所減小，說(shuō)明升力有所減小。這是由于結(jié)冰改變了槳葉外形，導(dǎo)致冰形后方出現(xiàn)了明顯的負(fù)壓區(qū)，使得氣動(dòng)特性惡化。從圖 6和圖 7中的壓力分布也可以看出結(jié)冰前槳葉表面的形狀尚未被破壞，此時(shí)流線(xiàn)緊貼壁面。結(jié)冰后由于槳葉表面的形狀被破壞，在前緣冰形后方出現(xiàn)一個(gè)負(fù)壓區(qū)域，流線(xiàn)中也可以看出此處存在渦，這導(dǎo)致槳葉阻力增加，升力減小，且這種趨勢(shì)隨著展向位置逐漸增強(qiáng)。

圖5 不同截面處壓力系數(shù)曲線(xiàn)對(duì)比

圖7 結(jié)冰后各個(gè)截面壓力分布情況

5 結(jié)語(yǔ)

本文對(duì)直升機(jī)旋翼槳葉結(jié)冰后氣動(dòng)特性進(jìn)行了研究，首先給出了一套三維旋翼槳葉結(jié)冰數(shù)值模擬方法，對(duì)其中的流場(chǎng)計(jì)算、水滴運(yùn)動(dòng)、結(jié)冰生成和幾何重構(gòu)進(jìn)行介紹。隨后計(jì)算了結(jié)冰前后旋翼槳葉的氣動(dòng)特性，分析了結(jié)冰對(duì)三維槳葉氣動(dòng)性能的影響。計(jì)算結(jié)果表明：

1)本文方法可以針對(duì)三維旋翼槳葉定量分析結(jié)冰對(duì)氣動(dòng)特性的影響。結(jié)冰會(huì)使槳葉整體升力下降，阻力上升，懸停效率下降超過(guò)20%。

2)結(jié)冰后旋翼槳葉各個(gè)截面的升力系數(shù)均發(fā)生降低，前緣處壓力系數(shù)波動(dòng)比較大，壓力系數(shù)會(huì)有一個(gè)突變的過(guò)程，阻力明顯增加。

3)結(jié)冰后槳葉原有氣動(dòng)外形被破壞，氣流流線(xiàn)發(fā)生分離，冰形后方甚至產(chǎn)生了渦，導(dǎo)致了氣動(dòng)特性的惡化。