民用飛機(jī)機(jī)頭冰脫落特性數(shù)值模擬

馮麗娟，譚廣琨，張 強(qiáng)，易 賢

(1.上海飛機(jī)設(shè)計(jì)研究院總體氣動(dòng)部，上海 201210; 2.中國(guó)空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川綿陽(yáng) 621000)

民用飛機(jī)機(jī)頭冰脫落特性數(shù)值模擬

馮麗娟1，*，譚廣琨1，張 強(qiáng)1，易 賢2

(1.上海飛機(jī)設(shè)計(jì)研究院總體氣動(dòng)部，上海 201210; 2.中國(guó)空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川綿陽(yáng) 621000)

建立了帶動(dòng)力條件下飛機(jī)表面冰塊脫落的數(shù)值計(jì)算方法，包括全機(jī)帶動(dòng)力的空氣流場(chǎng)計(jì)算方法、冰塊運(yùn)動(dòng)中氣動(dòng)力和力矩的確定方法以及冰塊運(yùn)動(dòng)的六自由度方程數(shù)值求解方法。采用以上方法對(duì)某民用支線飛機(jī)巡航和進(jìn)場(chǎng)狀態(tài)下的機(jī)頭冰脫落特性進(jìn)行了模擬，對(duì)比了迎角對(duì)冰脫落的影響，分析了冰脫落的速度和軌跡關(guān)系，得到了冰塊被吸入發(fā)動(dòng)機(jī)的概率，為發(fā)動(dòng)機(jī)吞冰設(shè)計(jì)、飛機(jī)冰脫落的適航符合性驗(yàn)證以及自然結(jié)冰試飛提供了較好的支持。

民用飛機(jī);機(jī)頭結(jié)冰;冰脫落;計(jì)算流體力學(xué);數(shù)值模擬

0 引言

飛機(jī)結(jié)冰廣泛存在于飛行過(guò)程中，是造成危害飛行安全的主要隱患之一［1-3］。飛機(jī)結(jié)冰主要發(fā)生在機(jī)翼和尾翼前緣、發(fā)動(dòng)機(jī)唇口、風(fēng)擋和機(jī)頭等部位［4］。除了機(jī)頭，以上部位通常會(huì)采取防除冰措施［5］。未采取防護(hù)措施的機(jī)頭結(jié)冰后，如果發(fā)生冰脫落現(xiàn)象，會(huì)對(duì)飛行安全產(chǎn)生影響。例如，脫落冰塊打到后機(jī)身或尾翼，可能損壞飛機(jī)部件，尤其當(dāng)脫落冰塊的質(zhì)量大于發(fā)動(dòng)機(jī)吞冰能力限制時(shí)，被發(fā)動(dòng)機(jī)吸入后會(huì)打壞發(fā)動(dòng)機(jī)葉片，造成發(fā)動(dòng)機(jī)嚴(yán)重推力損失甚至停機(jī)等后果。如何表明在各種飛行條件下，機(jī)頭上的結(jié)冰脫落后對(duì)飛機(jī)發(fā)動(dòng)機(jī)以及機(jī)體的安全不造成影響，滿足相關(guān)適航條款要求，對(duì)于該民機(jī)的適航取證是一個(gè)難點(diǎn)。如果在風(fēng)洞中進(jìn)行冰塊投放試驗(yàn)，一方面試驗(yàn)費(fèi)用較高，無(wú)法覆蓋大量的飛行條件，同時(shí)，現(xiàn)階段還不能做到完全真實(shí)模擬，尤其是無(wú)法在試驗(yàn)中評(píng)估發(fā)動(dòng)機(jī)動(dòng)力對(duì)冰脫落特性的影響。因此，有必要發(fā)展冰脫落研究的數(shù)值方法。

國(guó)外對(duì)于冰脫落研究的公開(kāi)報(bào)道不多，主要是用數(shù)值手段［6-9］或?qū)嶒?yàn)［10］分析冰塊的運(yùn)動(dòng)軌跡，研究通常采用簡(jiǎn)化的方塊冰形，大部分研究未考慮飛機(jī)帶動(dòng)力的情況。國(guó)內(nèi)近年來(lái)才有少量的研究。文獻(xiàn)［11］研究了冰形外部壓力和冰層與飛機(jī)蒙皮之間剪切應(yīng)力對(duì)冰脫落的影響，探討了冰脫落發(fā)生的原因和條件;文獻(xiàn)［12］采用四自由度模型計(jì)算分析了不同尺寸、形狀冰粒的脫落軌跡，未考慮三維流場(chǎng)和冰塊運(yùn)動(dòng)的六自由度;文獻(xiàn)［13］計(jì)算模擬了全機(jī)流場(chǎng)中冰塊的六自由度運(yùn)動(dòng)，但由于沒(méi)有考慮到發(fā)動(dòng)機(jī)的動(dòng)力影響，對(duì)近發(fā)動(dòng)機(jī)區(qū)域脫落冰塊的運(yùn)動(dòng)軌跡的分析與真實(shí)情況下脫落冰塊的軌跡會(huì)存在較大差異。

本文建立了一種考慮飛機(jī)發(fā)動(dòng)力影響的、基于六自由度運(yùn)動(dòng)方程的飛機(jī)冰脫落數(shù)值模擬方法，并對(duì)采用尾吊式發(fā)動(dòng)機(jī)布局的某民用支線飛機(jī)的機(jī)頭冰脫落特性進(jìn)行了計(jì)算模擬，得到了冰塊運(yùn)動(dòng)軌跡以及冰塊被吸入發(fā)動(dòng)機(jī)的概率。

1 冰脫落計(jì)算思路［13］

脫落冰塊的運(yùn)動(dòng)是一個(gè)非定常過(guò)程，同時(shí)與飛機(jī)流場(chǎng)存在相互干擾作用，考慮到冰塊尺寸遠(yuǎn)小于飛機(jī)，因此忽略冰塊對(duì)飛機(jī)流場(chǎng)的干擾效果，只考慮飛機(jī)流場(chǎng)對(duì)冰塊的作用，計(jì)算思路為:

(1)選取典型外形的冰塊;

(2)計(jì)算冰塊在各種姿態(tài)下的氣動(dòng)特性，建立冰塊的氣動(dòng)數(shù)據(jù)庫(kù);

(3)采用CFD方法計(jì)算飛機(jī)三維流場(chǎng)，獲得飛機(jī)流場(chǎng)分布;

(4)在上述工作基礎(chǔ)上，求解剛體運(yùn)動(dòng)六自由度方程，獲得脫落冰塊的運(yùn)動(dòng)軌跡;

(5)根據(jù)軌跡計(jì)算結(jié)果，分析冰塊被吸入發(fā)動(dòng)機(jī)的概率。

2 流場(chǎng)模擬方法

通過(guò)求解附帶k-ω SST湍流模型的RANS方程獲得飛機(jī)繞流流場(chǎng)以及冰塊隨迎角和側(cè)滑角變化的氣動(dòng)力數(shù)據(jù)。無(wú)量綱參數(shù)采用:密度、速度、壓力、能量、溫度 γ ( γ－1)、長(zhǎng)度 L、時(shí)間L/和粘性系數(shù)。三維直角坐標(biāo)系下守恒型RANS方程可以表達(dá)為:

采用LU-SGS格式的隱式時(shí)間推進(jìn)法求解控制方程，對(duì)流項(xiàng)離散采用Jameson中心格式，粘性項(xiàng)采用中心離散［14］。

飛機(jī)流場(chǎng)計(jì)算考慮了帶動(dòng)力特性的復(fù)雜構(gòu)型模擬，這樣能夠更為準(zhǔn)確地模擬飛機(jī)實(shí)際的流場(chǎng)特征，使得發(fā)動(dòng)機(jī)前方的冰脫落軌跡更加符合實(shí)際情況，提高了冰脫落軌跡計(jì)算的準(zhǔn)確性。

典型的渦輪風(fēng)扇發(fā)動(dòng)機(jī)計(jì)算模型如圖1所示，對(duì)發(fā)動(dòng)機(jī)進(jìn)排氣性能利用特定的計(jì)算邊界條件加以控制，其中風(fēng)扇入口為進(jìn)氣口，風(fēng)扇出口(外涵道)和主噴口(內(nèi)涵道)為排氣口。將進(jìn)氣口作為流場(chǎng)出口，給定質(zhì)量流量，按準(zhǔn)一維等熵流理論計(jì)算其他變量;排氣口作為流場(chǎng)入口，給定流動(dòng)方向、總溫和總壓。

圖1 渦輪風(fēng)扇發(fā)動(dòng)機(jī)計(jì)算模型示意圖Fig.1 Turbo fan engine simplified model

3 冰塊運(yùn)動(dòng)軌跡計(jì)算方法

冰塊運(yùn)動(dòng)軌跡通過(guò)數(shù)值求解剛體運(yùn)動(dòng)的六自由度方程獲得［15］。在六自由度方程求解中，以飛機(jī)的體軸系為參考來(lái)描述脫落物的質(zhì)心運(yùn)動(dòng)速度、質(zhì)心位移，冰塊脫落前的體軸系與飛機(jī)的體軸系之間存在俯仰角、偏航角和滾轉(zhuǎn)角三個(gè)初始?xì)W拉角。采用四階龍格-庫(kù)塔方法求解，具體步驟為:

(1)從脫落瞬間開(kāi)始，在軌跡模擬的任一時(shí)刻，利用冰塊的質(zhì)心坐標(biāo)從飛機(jī)流場(chǎng)數(shù)據(jù)中讀取所在位置的氣流速度，再根據(jù)冰塊姿態(tài)求得相對(duì)來(lái)流迎角和側(cè)滑角，進(jìn)而從冰塊氣動(dòng)數(shù)據(jù)庫(kù)中插值得到氣動(dòng)力和力矩;

(2)以該時(shí)刻的氣動(dòng)力和力矩作為輸入條件，求解剛體運(yùn)動(dòng)的六自由度方程，得到冰塊在下一時(shí)刻的位置和姿態(tài);

(3)根據(jù)冰塊新的位置和姿態(tài)，進(jìn)行下一時(shí)間步的計(jì)算，如此反復(fù)迭代，直到冰塊運(yùn)動(dòng)至特定區(qū)域。

4 民用飛機(jī)機(jī)頭冰脫落計(jì)算分析

4.1 脫落冰塊外形的選取

飛機(jī)機(jī)頭結(jié)冰外形通常呈圓盤(pán)狀，如圖2所示。

圖2 飛機(jī)機(jī)頭結(jié)冰外形Fig.2 Ice shape on nose

選取尺寸為1/3圓盤(pán)的扇形冰塊表征機(jī)頭圓盤(pán)冰破裂后形成的脫落冰塊，如圖3所示。

4.2 冰塊氣動(dòng)數(shù)據(jù)庫(kù)的建立

圖3 扇形冰塊Fig.3 Fan like ice shape

冰塊氣動(dòng)數(shù)據(jù)庫(kù)在冰塊軌跡計(jì)算過(guò)程中提供氣動(dòng)力和力矩信息，包含三個(gè)笛卡爾坐標(biāo)軸方向的氣動(dòng)力系數(shù)和力矩系數(shù)。圖4為脫落冰塊的計(jì)算網(wǎng)格，遠(yuǎn)場(chǎng)邊界為10倍冰塊長(zhǎng)度。

圖4 冰塊空間網(wǎng)格及遠(yuǎn)場(chǎng)邊界Fig.4 Mesh of ice and far field

圖5為冰塊氣動(dòng)力隨迎角和側(cè)滑角變化的數(shù)據(jù)庫(kù)，其中橫坐標(biāo)為迎角，變化區(qū)間為(-180°，180°)，不同的曲線族按側(cè)滑角分類(lèi)，變化區(qū)間為(-90°，90°)。

4.3 計(jì)算條件

以某尾吊式布局的民用支線飛機(jī)作為研究對(duì)象，選用的扇形脫落冰塊半徑為0.12 m，厚度0.04 m，冰塊的質(zhì)量為0.454 kg。

機(jī)頭脫落點(diǎn)位置:以機(jī)頭駐點(diǎn)為圓心、半徑為0.25 m的圓上選取均勻分布的6個(gè)點(diǎn)作為機(jī)頭冰塊的初始脫落位置，如圖6所示。

圖5 冰塊氣動(dòng)數(shù)據(jù)庫(kù)Fig.5 Ice slice aerodynamic database

圖6 冰塊在機(jī)頭初始脫落位置Fig.6 Initial position of ice slice on nose

冰塊初始姿態(tài)表示為不同的歐拉角(俯仰角、偏航角、滾轉(zhuǎn)角)的組合。

飛機(jī)的計(jì)算網(wǎng)格采用結(jié)構(gòu)網(wǎng)格，物面第一層為平均弦長(zhǎng)的10-5，遠(yuǎn)場(chǎng)邊界為10倍飛機(jī)長(zhǎng)度，網(wǎng)格布局注意機(jī)身的保形，同時(shí)在流場(chǎng)變化較快的地方進(jìn)行局部加密。

計(jì)算了巡航和進(jìn)場(chǎng)兩種狀態(tài)，具體計(jì)算狀態(tài)見(jiàn)表1。

表1 飛機(jī)飛行狀態(tài)Table 1 Flight condition

4.4 巡航狀態(tài)冰脫落特性計(jì)算分析

圖7為巡航狀態(tài)不同迎角下脫落冰塊軌跡的前側(cè)視圖，曲線顏色代表冰塊質(zhì)心的移動(dòng)速度，由于馬赫數(shù)較大，脫落冰塊受氣動(dòng)力影響較重力更明顯，冰塊群整體位于發(fā)動(dòng)機(jī)上側(cè)，冰塊運(yùn)動(dòng)到發(fā)動(dòng)機(jī)唇口截面時(shí)的最大速度達(dá)到38 m/s。圖8為飛機(jī)左側(cè)半模冰塊軌跡與發(fā)動(dòng)機(jī)唇口軸向截面的交點(diǎn)分布圖，由圖可知，隨著迎角的增大，軌跡交點(diǎn)整體上移，其中4.2°迎角時(shí)吸入發(fā)動(dòng)機(jī)的點(diǎn)數(shù)最多，概率為4.32%，隨著迎角增大，吸入發(fā)動(dòng)機(jī)的冰塊數(shù)明顯減少，迎角6.6°時(shí)，所計(jì)算的脫落軌跡中，只有一條進(jìn)入發(fā)動(dòng)機(jī)唇口，概率為0.62%，而迎角8.4°則未發(fā)現(xiàn)吸入發(fā)動(dòng)機(jī)的冰塊。

圖7 巡航構(gòu)型冰脫落軌跡Fig.7 Ice shedding trajectories of cruise case

4.5 進(jìn)場(chǎng)狀態(tài)冰脫落特性計(jì)算分析

圖9為進(jìn)場(chǎng)狀態(tài)不同迎角下脫落冰塊軌跡的前側(cè)視圖，與巡航狀態(tài)相比，由于馬赫數(shù)明顯減小，氣動(dòng)力的作用減弱，因此脫落軌跡向下偏移較多，冰塊運(yùn)動(dòng)到發(fā)動(dòng)機(jī)唇口截面時(shí)的最大速度為30 m/s。圖10為飛機(jī)左側(cè)半模冰塊軌跡與發(fā)動(dòng)機(jī)唇口軸向截面的交點(diǎn)分布圖，在4.3°迎角時(shí)軌跡交點(diǎn)均位于發(fā)動(dòng)機(jī)下部，沒(méi)有冰塊打入發(fā)動(dòng)機(jī)，隨著迎角增大，冰塊軌跡上揚(yáng)，但吸入發(fā)動(dòng)機(jī)的冰塊仍然很少，迎角為7.3°和10.3°時(shí)，分別只有一條軌跡打入發(fā)動(dòng)機(jī)，概率為0.62%。

圖8 巡航構(gòu)型冰脫落軌跡與唇口軸向截面的交點(diǎn)Fig.8 Footprint of shedding trajectories on engine plane of cruise case

圖9 進(jìn)場(chǎng)構(gòu)型冰脫落軌跡Fig.9 Ice shedding trajectories of approach case

圖10 巡航構(gòu)型冰脫落軌跡與唇口軸向截面的交點(diǎn)Fig.10 Footprint of shedding trajectories on engine plane of approach case

5 結(jié)論

本文建立了冰脫落計(jì)算分析方法，并針對(duì)某型民用支線飛機(jī)帶動(dòng)力條件下機(jī)頭冰脫落問(wèn)題進(jìn)行了計(jì)算分析，得到如下結(jié)論:

(1)冰塊運(yùn)動(dòng)軌跡受重力，來(lái)流速度和迎角的影響，正迎角下氣動(dòng)力使軌跡上揚(yáng)，重力則使軌跡下偏，隨著來(lái)流速度以及迎角的增大，氣動(dòng)力作用影響加大，冰塊軌跡升高;

(2)通過(guò)對(duì)該民機(jī)機(jī)頭冰脫落的軌跡被吸入發(fā)動(dòng)機(jī)的概率分析來(lái)看，帶動(dòng)力條件下，對(duì)于本文所計(jì)算的狀態(tài)，該尾吊式發(fā)動(dòng)機(jī)吸入機(jī)頭脫落冰的概率較低。

(3)本文的數(shù)值方法可以為民用飛機(jī)冰脫落適航符合性驗(yàn)證提供支持，也可作為評(píng)估自然結(jié)冰試飛中冰脫落對(duì)發(fā)動(dòng)機(jī)安全影響的手段。

［1］ Cebeci T，Kafyeke F.Aircraft icing［J］.Annu.Rev.Fluid Mech，2003，35:11-21.

［2］ Kind R J，Potapczuk M G.Experimental and computational simulation of in-flight icing phenomena［J］.Progress in Aerospace Science，1998，34:275-345.

［3］ Frank T L，Abdollah K.Effects of ice accretions on aircraft aer-odynamics［J］.Progress in Aerospace Sciences，2001，37(8):669-767.

［4］ Yi Xian，Zhu Guolin，et al.Numerically simulation of ice accretion on airfoil and research of icing wind tunnel test scale rule［J］.Acta Aerodynamica Sinica，2002，20(4):428-433.(in Chinese)

易賢，朱國(guó)林，等.飛機(jī)積冰的數(shù)值計(jì)算與積冰試驗(yàn)相似準(zhǔn)則研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2002，20(4):428-433.

［5］ Qiu Xiegang，Han Fenghua.Aircraft anti-icing system［M］.Beijing:Aviation Professional Teaching Material Compliance Group，1985.(in Chinese)

裘燮剛，韓鳳華.飛機(jī)防冰系統(tǒng)［M］.北京:航空專業(yè)教材編審組，1985.

［6］ Luis C de Castro Santos，Ramon Papa，Marcello do Areal S.Ferrari.A simulation model for ice impact risk evaluation［R］.AIAA 2003-0030.

［7］ Wright W B，Keith T G，DeWitt K J.Numerical simulation of icing，deicing and shedding［R］.AIAA-91-0665，1991.

［8］ Michael Papadakis，Hsiung-Wei Yeong，Ian G Suares，et al.Experimental and computational investigation of ice shedding from aircraft surfaces［R］.AIAA 2006-1010，2006.

［9］ Hsiung-Wei Yeong，Michael Papadakis，Koji Shimoi.Ice trajectories and Monte Carlo analyses for a business jet［R］.AIAA 2009-3973.

［10］ Michael Papadakis，Hsiung-Wei Yeong，Koji Shimoi，et al.Ice shedding experiments with simulated ice shapes［R］.AIAA 2009-3972.

［11］Xiao Chunhua，Gui Yewei，Yi Xian，et al.Pressure distribution and stress calculation of glaze ice on iced airfoil［J］.Journal of Aerospace Power，2009，24(7):1457-1463.(in Chinese)

肖春華，桂業(yè)偉，易賢，等.結(jié)冰翼型表面明冰的壓力分布和應(yīng)力計(jì)算［J］.航空動(dòng)力學(xué)報(bào)，2009，24(7):1457-1463.

［12］ Tan Meijing，Sun Zhiguo，Zhu Chunling，et al.Four degree of Freedom simulation of ice shedding for aircraft icing［J］.Aeronautical Computing Technique，2012，42(2):100-104.(in Chinese)

檀妹靜，孫志國(guó)，朱春玲.四自由度飛機(jī)結(jié)冰脫落模擬研究［J］.航空計(jì)算技術(shù)，2012，42(2):100-104.

［13］Wang Jiantao，Yi Xian，Xiao Zhongyun，et al.Numerical simulation of ice shedding from ARJ21-700［J］.Acta Aerodynamica Sinica，2013，31(4):430-436.(in Chinese)

王建濤，易賢，肖中云，等.ARJ21-700飛機(jī)冰脫落數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2013，31(4):430-436.

［14］ Mou Bin，Xiao Zhongyun，Liu Gang，et al.Low-speed preconditioning based on Roe scheme［J］.Acta Aerodynamica Sinica，2013，31(4):430-436.(in Chinese)

牟斌，肖中云，劉剛，等.基于Roe格式的低速預(yù)處理方法研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2010，28(2):125-131.

［15］Xu Mingyou，Ding Songbin.Flight dynamic［M］.Beijing:Science Press，2003.(in Chinese)

徐明友，丁松濱.飛行動(dòng)力學(xué)［M］.北京:科學(xué)出版社，2003.

Numerical simulation of ice shedding from civil airplane nose

Feng Lijuan1，*，Tan Guangkun1，Zhang Qiang1，Yi Xian2
(1.General and Aerodynamics Department，Shanghai Aircraft Design Research Institute，Shanghai 201210，China; 2.State Key Laboratory of Aerodynamics，China Aerodynamics Research and Development Center，Mianyang 621000，China)

The shedding ice blocks from a civil airplane nose seriously threat the safety of the civil airplane，and as result it is very important to thoroughly investigate the shedding characteristics.A numerical simulation method of ice shedding was proposed to simulate the trajectories of the shedding ice blocks from airplane nose by establishing the force and moment database of the shedding ice blocks using CFD(Computational Fluid Dynamics)tool，and solving the six degree-of-freedom equations numerically，the effect of engine power was also considered.Employing this method，the trajectories of the shedding ice from a civil airplane nose with aft-mounted engines were simulated，the effect of angle of attack of airplane on the trajectories of shedding ice were studied，the relationship among the velocity，the trajectories of the shedding ice and the probability of shedding ice blocks ingested into the engines was analyzed.The research work is of help for the design of engine ingestion ability，the compliance work of nose ice shedding and the nature icing flight test of the civil aircraft.

civil airplane;airplane nose icing;ice shedding;computational fluid dynamics;numerical simulation

V244.1

A

10.7638/kqdlxxb-2014.0005

0258-1825(2015)04-0536-06

2014-02-17;

2014-05-14

國(guó)家自然科學(xué)基金(11172314)

馮麗娟*(1976-)，女，碩士，高工，主要從事民機(jī)結(jié)冰適航取證方面的研究工作.E-mail:fenglijuan@comac.cc

馮麗娟，譚廣琨，張強(qiáng)，等.民用飛機(jī)機(jī)頭冰脫落特性數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2015，33(4):536-541.

10.7638/kqdlxxb-2014.0005 Feng L J，Tan G K，Zhang Q，et al.Numerical simulation of ice shedding from civil airplane nose［J］.Acta Aerodynamica Sinica，2015，33(4):536-541.