直升機旋翼結(jié)冰后三維槳葉氣動特性數(shù)值分析

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

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

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

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

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

1 計算模型

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

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

圖1 旋翼槳葉流場計算域

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

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

2.1 空氣流場計算

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

(1)

(2)

(3)

2.2 水滴撞擊特性計算

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

(4)

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

(5)

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

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

(6)

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

2.3 結(jié)冰生成計算

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

(7)

(8)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 不同截面處壓力系數(shù)曲線對比

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

5 結(jié)語

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

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

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

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