常士楠，楊秋明，李 延

(北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，北京 100191)

0 引言

飛機(jī)在結(jié)冰氣象條件下飛行時(shí)會(huì)發(fā)生結(jié)冰現(xiàn)象。飛機(jī)機(jī)翼表面前緣形成的冰層會(huì)改變機(jī)翼的氣動(dòng)外形，惡化機(jī)翼的氣動(dòng)特性，導(dǎo)致升力減小、阻力增大、失速迎角減小、力矩特性變差，嚴(yán)重危及整架飛機(jī)的飛行安全。因此對(duì)機(jī)翼的結(jié)冰問(wèn)題進(jìn)行深入研究非常必要。歐美和前蘇聯(lián)的相關(guān)研究機(jī)構(gòu)對(duì)飛機(jī)結(jié)冰問(wèn)題進(jìn)行了長(zhǎng)期而大量的細(xì)致研究，在試驗(yàn)研究和數(shù)值模擬方面都取得了豐碩成果。特別是近年來(lái)隨著計(jì)算流體力學(xué)(CFD)技術(shù)的發(fā)展，通過(guò)數(shù)值模擬手段對(duì)結(jié)冰過(guò)程和結(jié)冰后氣動(dòng)特性進(jìn)行數(shù)值預(yù)測(cè)成為廣大研究者關(guān)注的熱點(diǎn)問(wèn)題［1-2］。雖然國(guó)內(nèi)對(duì)機(jī)翼結(jié)冰問(wèn)題的研究起步較晚，但是近年來(lái)越來(lái)越多的學(xué)者開(kāi)始從事這方面的工作，也取得了不少成果:南航張大林提出一種邊界移動(dòng)技術(shù)來(lái)模擬翼型表面霜冰的生長(zhǎng)過(guò)程［3］，西工大李鳳蔚研究了時(shí)間步長(zhǎng)和網(wǎng)格數(shù)對(duì)冰形的影響［4］，空氣動(dòng)力研究與發(fā)展中心的朱國(guó)林和易賢對(duì)二維、三維外形的機(jī)翼前緣積冰進(jìn)行了數(shù)值模擬，同時(shí)對(duì)積冰試驗(yàn)相似準(zhǔn)則進(jìn)行了研究［5］。

對(duì)冰形進(jìn)行準(zhǔn)確預(yù)測(cè)的難點(diǎn)在于翼型表面產(chǎn)生一定厚度的冰層后壁面邊界條件發(fā)生變化，網(wǎng)格、外流場(chǎng)、水滴撞擊特性都將隨之改變。本文利用Fluent流體計(jì)算商業(yè)軟件中的動(dòng)網(wǎng)格技術(shù)對(duì)翼型表面結(jié)冰的整個(gè)過(guò)程進(jìn)行了模擬，在每個(gè)時(shí)間步長(zhǎng)內(nèi)完成網(wǎng)格隨著壁面邊界的移動(dòng)而更新、周?chē)鲌?chǎng)和水滴撞擊特性的重新計(jì)算、結(jié)冰冰形計(jì)算及壁面邊界的重構(gòu)工作，如此循環(huán)直至所需的結(jié)冰時(shí)間，得到最終的冰形。本文的計(jì)算方法和程序適用于對(duì)各種冰形的模擬。

1 流場(chǎng)計(jì)算

1.1 網(wǎng)格劃分

網(wǎng)格生成是對(duì)流場(chǎng)進(jìn)行數(shù)值計(jì)算的基礎(chǔ)，本文以NACA0012翼型為例，根據(jù)Fluent軟件中動(dòng)網(wǎng)格功能的要求和結(jié)冰實(shí)際情況采用非結(jié)構(gòu)化的三角形網(wǎng)格，便于對(duì)網(wǎng)格的局部重構(gòu)和光順，翼型壁面附近特別是前緣部分網(wǎng)格密度增大用來(lái)捕捉附面層區(qū)域的流動(dòng)信息，網(wǎng)格節(jié)點(diǎn)數(shù)為111816，網(wǎng)格單元數(shù)為188738，網(wǎng)格劃分如圖1所示:

圖1 NACA0012翼型網(wǎng)格示意圖Fig.1 Diagram of grid around NACA0012 airfoil

1.2 流場(chǎng)計(jì)算

根據(jù)質(zhì)量、動(dòng)量、能量三大守恒關(guān)系建立二維非定常Navier-Stokes方程，針對(duì)某一控制體Ω寫(xiě)成如下積分形式［4］:

式中V為微元控制體，U為質(zhì)量、動(dòng)量和能量各守恒變量，s為翼型的壁面邊界，F(xiàn)為無(wú)粘通量，F(xiàn)v為粘性通量，n為控制體的外法線向量，Re為翼型的特征雷諾數(shù)。

1.3 水滴撞擊特性

局部收集系數(shù)是結(jié)冰模擬過(guò)程中的一個(gè)重要參數(shù)，本文采用拉格朗日軌跡追蹤法通過(guò)求解水滴運(yùn)動(dòng)方程來(lái)跟蹤各個(gè)水滴的運(yùn)動(dòng)軌跡，進(jìn)而確定機(jī)翼壁面上局部收集系數(shù)的分布。由于空氣中液態(tài)水含量在克量級(jí)，水滴的平均容積直徑在微米量級(jí)，因此可作如下合理假設(shè)［6］:(1)水滴在運(yùn)動(dòng)過(guò)程中和空氣不發(fā)生熱交換、不蒸發(fā)，水滴的物性參數(shù)不變化;(2)水滴在運(yùn)動(dòng)過(guò)程中不凝結(jié)、不破碎、不變形;(3)空氣中水滴的存在不影響空氣的流動(dòng)，空氣的紊流和脈動(dòng)不影響水滴的運(yùn)動(dòng)，水滴的初速度與自由來(lái)流的速度相等。

作用在水滴上的力僅有粘性阻力、重力和空氣浮力，根據(jù)牛頓第二定律，水滴軌跡運(yùn)動(dòng)方程可寫(xiě)為:

水滴的運(yùn)動(dòng)方程可改寫(xiě)為:

式中，d為水滴直徑。

將前面流場(chǎng)速度分布的計(jì)算結(jié)果作為初始條件，上式的求解可看成一個(gè)一階常微分方程的初值問(wèn)題，因此采用四步Runge-Kutta法求解。得到各個(gè)水滴的軌跡后就可確定翼型壁面各控制體上的局部收集系數(shù)。

從光潔翼型開(kāi)始到每隔單個(gè)時(shí)間步長(zhǎng)翼型周?chē)牧鲌?chǎng)計(jì)算結(jié)束后，都要采用上述方法計(jì)算局部收集系數(shù)的分布，這樣才能體現(xiàn)出結(jié)冰后壁面邊界變化對(duì)局部收集系數(shù)分布的影響。

1.4 對(duì)流換熱系數(shù)的計(jì)算

翼型表面的對(duì)流換熱系數(shù)對(duì)最終的冰形具有重要影響。傳統(tǒng)的計(jì)算對(duì)流換熱系數(shù)的方法有雷諾比擬［8］和光表面的邊界層積分法［5］。雖然這兩種方法使用簡(jiǎn)單，但沒(méi)有考慮結(jié)冰后翼型表面粗糙度對(duì)對(duì)流換熱系數(shù)的影響，與實(shí)際情況有明顯差異。本文采用Kays和Crawford提出的考慮粗糙度的積分邊界層法計(jì)算對(duì)流換熱系數(shù)［9］，與實(shí)際的結(jié)冰情況更加接近。

目前對(duì)結(jié)冰翼型表面粗糙度形成的機(jī)理還不清楚，粗糙度分布又非常復(fù)雜，因此很難進(jìn)行精確模擬。本文采用等效的沙粒粗糙度來(lái)考慮粗糙度對(duì)積冰過(guò)程特別是對(duì)流換熱系數(shù)的影響［10］。等效沙粒粗糙度高度ks的計(jì)算公式為:

其中LWC為液態(tài)水含量(g/m3)，T∞為環(huán)境溫度，V∞為飛行速度，c為翼型的弦長(zhǎng)。

1.4.1 層流區(qū)域

在層流區(qū)域，局部位置的斯坦頓數(shù)［3］可用下式表示:

其中 Ge為邊界層邊緣的質(zhì)量流速，Ge=ρUe;C1、C2、C3是與 Pr相關(guān)的函數(shù)，Pr取常數(shù)0.72，C1取0.41，C2取0.44，C3取 1.88［11］。根據(jù) Chilton-Colburn［9］分析，結(jié)冰表面對(duì)流換熱系數(shù)為:

式中ρa(bǔ)為空氣密度，cpa為空氣比熱，ue為邊界層邊緣速度，定義為99%的來(lái)流速度，St為無(wú)量綱斯坦頓數(shù)。

根據(jù)上式可推得層流區(qū)域的局部對(duì)流換熱系數(shù):

式中λ為空氣的導(dǎo)熱系數(shù)，υ為空氣的運(yùn)動(dòng)粘度。

1.4.2 層流到湍流過(guò)渡位置的確定

本文在考慮粗糙度因素的影響時(shí)，采用Von Doen hoff和Horton［12］的標(biāo)準(zhǔn)，即滿(mǎn)足粗糙度臨界雷諾數(shù)

時(shí)，流動(dòng)邊界層處于湍流區(qū)域。式中uk為y=ks處的速度，ks為等效沙粒粗糙高度［8］。而根據(jù)文獻(xiàn)［9］可知:

1.4.3 湍流區(qū)域

同層流區(qū)域情況類(lèi)似，首先求出邊界層局部位置的斯坦頓數(shù)，是等效砂粒粗糙度、剪力速度、摩擦系數(shù)和空氣參數(shù)的函數(shù)。即:

其中，湍流普朗特?cái)?shù)Prt取常數(shù)值0.9，Pr為0.72。式中Stk為無(wú)量綱粗糙度斯坦頓數(shù):

式中C1取常數(shù)值1.0［14］，u為剪力速度，不同文獻(xiàn)采用不同表達(dá)式，本文采用Kays及 Crawford［9］中表達(dá)式:

式中Cf為邊界層區(qū)域局部位置的摩擦系數(shù)，Cf的表達(dá)式為［5］:

θturb為湍流邊界層的動(dòng)量厚度:

式中str是指壁面上駐點(diǎn)到轉(zhuǎn)捩點(diǎn)的弧長(zhǎng)，θstr是指轉(zhuǎn)捩點(diǎn)處層流邊界層的動(dòng)量厚度。

由此可得到無(wú)量綱粗糙度斯坦頓數(shù)Stk:

將上式帶入St的表達(dá)式并簡(jiǎn)化，得:

則邊界層湍流區(qū)域的局部對(duì)流換熱系數(shù)的表達(dá)式可表示為:

每次流場(chǎng)重新計(jì)算后都要用上述方法重新計(jì)算對(duì)流換熱系數(shù)的分布，為下一步結(jié)冰熱力學(xué)模型的求解奠定基礎(chǔ)。

完成每個(gè)時(shí)間步長(zhǎng)的結(jié)冰計(jì)算后都需要在壁面邊界移動(dòng)后新生成的計(jì)算域內(nèi)重復(fù)上述流場(chǎng)計(jì)算、水滴撞擊特性計(jì)算和對(duì)流換熱系數(shù)計(jì)算。利用重新計(jì)算得到的上述結(jié)果進(jìn)行下一個(gè)時(shí)間步長(zhǎng)內(nèi)的結(jié)冰計(jì)算直到完成整個(gè)結(jié)冰時(shí)間內(nèi)的結(jié)冰過(guò)程的模擬。

2 結(jié)冰過(guò)程的數(shù)值模擬

2.1 Messinger積冰熱力學(xué)模型

目前對(duì)結(jié)冰問(wèn)題的研究大都采用經(jīng)典的Messinger結(jié)冰熱力學(xué)模型［15］，其基本思想是:選取結(jié)冰前機(jī)翼表面和已結(jié)冰后積冰表面的微元體為控制體，并假設(shè)為準(zhǔn)定常過(guò)程來(lái)建立質(zhì)量守恒和能量守恒方程，質(zhì)量守恒方程中包含撞擊水滴的質(zhì)量流、流入控制體水的質(zhì)量流、表面蒸發(fā)水的質(zhì)量流、凍結(jié)為冰的質(zhì)量流和流出控制體的質(zhì)量流，能量守恒方程中包含撞擊到表面的過(guò)冷水滴的能量mimiw，T、流入控制體的水滴的能量 miniw，sur(i-1)、流出控制體的水滴的能量moutiw，sur、蒸發(fā)水的能量 mevaiv，sur、凍結(jié)水的能量 mfizii，sur、氣流與積冰表面對(duì)流換熱釋放的能量qcΔs和控制體底部與壁面的導(dǎo)熱能量qcΔs，Δs為控制體的換熱面積。由質(zhì)量守恒和能量守恒關(guān)系分別建立方程:

控制體的質(zhì)量守恒和能量守恒關(guān)系具體見(jiàn)參考文獻(xiàn)［16］。

2.2 結(jié)冰熱力學(xué)模型的求解

為了便于求解Messinger結(jié)冰熱力學(xué)模型，引入凍結(jié)系數(shù)f，它表示每個(gè)控制體中結(jié)冰水的質(zhì)量與控制體總共收集水質(zhì)量的比值:

聯(lián)立質(zhì)量方程和能量方程，并將f帶入可得:

各物理參數(shù)的具體意義詳見(jiàn)文獻(xiàn)［17］，將各參數(shù)的表達(dá)式帶入式(23)后最終化簡(jiǎn)為含有表面溫度ts(℃)和凍結(jié)系數(shù)f的不定方程。需要增加由結(jié)冰過(guò)程中實(shí)際物理意義的約束條件進(jìn)行迭代求解。

2.3 壁面邊界重構(gòu)與光順

得到各控制體的結(jié)冰質(zhì)量mfrz后，可計(jì)算出每個(gè)時(shí)間步長(zhǎng)內(nèi)各控制體冰層的生長(zhǎng)高度Δh:

根據(jù)冰層始終沿壁面外法線方向生長(zhǎng)的假定，對(duì)翼型結(jié)冰后的壁面網(wǎng)格進(jìn)行重構(gòu)，

根據(jù)各控制體冰層生長(zhǎng)的高度和方向，將壁面上結(jié)冰區(qū)域的網(wǎng)格節(jié)點(diǎn)進(jìn)行移動(dòng)，設(shè)定網(wǎng)格單元的最小面積Smin，當(dāng)與移動(dòng)節(jié)點(diǎn)相關(guān)的網(wǎng)格單元面積小于Smin時(shí)就與相鄰網(wǎng)格單元合并，這樣可以保證網(wǎng)格移動(dòng)過(guò)程中三角形網(wǎng)格的畸變率不會(huì)增大，網(wǎng)格的質(zhì)量不會(huì)下降。這樣處理，只需要將壁面上結(jié)冰區(qū)域的節(jié)點(diǎn)和網(wǎng)格單元進(jìn)行移動(dòng)和重構(gòu)，翼型后緣無(wú)冰生成的區(qū)域和遠(yuǎn)場(chǎng)區(qū)域的網(wǎng)格不需要變化，因此可以避免對(duì)單個(gè)結(jié)冰時(shí)間結(jié)束后翼型邊界發(fā)生變化后的整個(gè)流場(chǎng)計(jì)算區(qū)域重新生成網(wǎng)格，提高效率。

冰層沿壁面外法線方向生長(zhǎng)可以得到結(jié)冰時(shí)間步長(zhǎng)后壁面各節(jié)點(diǎn)的新坐標(biāo)，但壁面上曲率變化劇烈的區(qū)域就會(huì)導(dǎo)致移動(dòng)后的節(jié)點(diǎn)光滑性變差，甚至出現(xiàn)相互交錯(cuò)的情況，嚴(yán)重影響網(wǎng)格質(zhì)量，因此在進(jìn)行動(dòng)網(wǎng)格過(guò)程之前需要在保持總結(jié)冰量不變的基礎(chǔ)上對(duì)壁面上移動(dòng)后的節(jié)點(diǎn)進(jìn)行光順。從駐點(diǎn)開(kāi)始，分別遍歷上下壁面上的節(jié)點(diǎn)，三個(gè)節(jié)點(diǎn)編號(hào)為i-1、i、i+1，每3個(gè)節(jié)點(diǎn)總能組成一個(gè)三角形或者形成一條直線，以i為頂點(diǎn)的夾角大于設(shè)定值時(shí)，則需將節(jié)點(diǎn)i進(jìn)行光順。光順后節(jié)點(diǎn)i的新坐標(biāo)為:

對(duì)各節(jié)點(diǎn)的光順需要遍歷多次，確保每個(gè)節(jié)點(diǎn)至少有一次作為頂點(diǎn)被光順過(guò)，同時(shí)要將駐點(diǎn)節(jié)點(diǎn)作為頂點(diǎn)進(jìn)行光順，確保上下壁面的光滑連接。根據(jù)網(wǎng)格節(jié)點(diǎn)的疏密程度，本文中夾角的設(shè)定值取150°。圖2是結(jié)冰時(shí)間為360s時(shí)經(jīng)過(guò)重構(gòu)與光順后的網(wǎng)格示意圖。

通過(guò)以上所述方法對(duì)壁面節(jié)點(diǎn)、網(wǎng)格的重構(gòu)和光順，可以保證冰形輪廓不失真，提高網(wǎng)格質(zhì)量，同時(shí)避免了對(duì)結(jié)冰后翼型邊界發(fā)生變化后的流場(chǎng)計(jì)算區(qū)域重新生成網(wǎng)格。完成網(wǎng)格的重構(gòu)和光順后，進(jìn)行下一個(gè)時(shí)間步長(zhǎng)的流場(chǎng)計(jì)算、水滴撞擊特性計(jì)算和結(jié)冰熱力學(xué)模型的求解，如此循環(huán)直至所需要的總的結(jié)冰時(shí)間。

3 算例分析

為了便于與冰風(fēng)洞實(shí)驗(yàn)結(jié)果比較，本文以NACA 0012翼型為例，分別模擬了環(huán)境溫度為 -6.1℃、-28.4℃和-13.4℃的三種典型冰形，計(jì)算條件如下:

表1 飛行條件和氣象參數(shù)表Table1 Flight condition and air parameter

圖3、圖4和圖5分別表示環(huán)境靜溫為-6.1℃時(shí)，結(jié)冰時(shí)間步長(zhǎng)為60s時(shí)翼型壁面上的局部收集系數(shù)、對(duì)流換熱系數(shù)和凍結(jié)系數(shù)的分布情況。由圖3、圖4和圖5可以看出，經(jīng)過(guò)每個(gè)結(jié)冰時(shí)間步長(zhǎng)后，翼型原有的外形發(fā)生變化，壁面上的局部收集系數(shù)和對(duì)流換熱系數(shù)隨之變化，特別是駐點(diǎn)附近的變化最為明顯。圖6表示了從潔凈翼型開(kāi)始每隔60s的翼型表面冰形的變化過(guò)程。環(huán)境靜溫較高時(shí)，駐點(diǎn)區(qū)域凍結(jié)系數(shù)較小，結(jié)冰量小，未凍結(jié)的液態(tài)水在氣流的吹拂下向下游溢流，最終形成槽狀冰。駐點(diǎn)區(qū)域冰形變化復(fù)雜，導(dǎo)致局部收集系數(shù)和對(duì)流換熱系數(shù)變化劇烈，這與圖3和圖4的結(jié)果也是相一致的。圖7表示了本文預(yù)測(cè)的最終冰形與LEWICE冰風(fēng)洞實(shí)驗(yàn)結(jié)果［16］和LEWICE軟件計(jì)算結(jié)果［17］的對(duì)比情況?？梢钥闯?，本文結(jié)果與實(shí)驗(yàn)結(jié)果基本一致，上下壁面處結(jié)冰極限基本與實(shí)驗(yàn)結(jié)果相吻合，只是上壁面處冰角的厚度大于實(shí)驗(yàn)結(jié)果，但冰形的總體效果好于LEWICE軟件計(jì)算結(jié)果。

3.2 環(huán)境溫度為-28.4℃時(shí)的楔形冰計(jì)算

圖8表示了環(huán)境靜溫為-28.4℃時(shí)經(jīng)過(guò)360s生成的楔形冰，可以看出本文預(yù)測(cè)的冰形與實(shí)驗(yàn)結(jié)果［16］吻合的很好，只是在下壁面處結(jié)冰極限小于實(shí)驗(yàn)結(jié)果，這可能是結(jié)冰熱力學(xué)模型中只考慮質(zhì)量和能量守恒關(guān)系，而沒(méi)有包含凍結(jié)過(guò)程中液態(tài)水在氣動(dòng)力作用下的運(yùn)動(dòng)規(guī)律。

圖3 不同結(jié)冰時(shí)間局部收集系數(shù)分布圖Fig.3 Local collection coefficient distribution of different icing time

圖4 不同結(jié)冰時(shí)間壁面對(duì)流換熱系數(shù)分布圖Fig.4 Wall convective heat-transfer coefficient distribution of different icing time

圖5 不同結(jié)冰時(shí)間凍結(jié)系數(shù)分布圖Fig.5 Freeze fraction distribution of different icing time

3.3 環(huán)境溫度為-13.4℃時(shí)的混合冰計(jì)算

圖9表示了環(huán)境靜溫為-13.4℃時(shí)經(jīng)過(guò)360s生成的典型的混合冰，可以看出本文結(jié)果在結(jié)冰極限位置與實(shí)驗(yàn)結(jié)果［16］基本一致，但是上壁面冰角處和下壁面溢流范圍處的冰形均有明顯差異，上壁面冰角的厚度小于實(shí)驗(yàn)結(jié)果?；旌媳嬗胁蹱畋托ㄐ伪奶攸c(diǎn)，預(yù)測(cè)結(jié)果與實(shí)驗(yàn)結(jié)果有較明顯差異，因此在結(jié)冰熱力學(xué)模型方面需要進(jìn)一步改進(jìn)，對(duì)結(jié)冰過(guò)程的機(jī)理和細(xì)節(jié)需要深入研究。

4 結(jié)論

圖6 不同結(jié)冰時(shí)間冰形圖Fig.6 Ice shapes of different icing time

圖7 環(huán)境溫度為-6.1℃時(shí)的冰形比較Fig.7 Comparison of ice shapes at temperature-6.1℃

圖8 環(huán)境溫度為-28.4℃時(shí)的冰形比較Fig.8 Comparison of ice shapes at temperature-28.4℃

圖9 環(huán)境溫度為-13.4℃時(shí)的冰形比較Fig.9 Comparison of ice shapes at temperature-13.4℃

利用Fluent動(dòng)網(wǎng)格模塊對(duì)翼型表面的結(jié)冰過(guò)程進(jìn)行了準(zhǔn)定常模擬，在每個(gè)結(jié)冰時(shí)間步長(zhǎng)內(nèi)完成網(wǎng)格隨著壁面邊界的移動(dòng)而更新、周?chē)鲌?chǎng)和水滴撞擊特性重新計(jì)算、對(duì)流換熱系數(shù)的重新計(jì)算結(jié)冰冰形計(jì)算及壁面邊界的重構(gòu)工作，最后計(jì)算了其它結(jié)冰條件相同而環(huán)境溫度不同的三種典型冰形，將預(yù)測(cè)結(jié)果與冰風(fēng)洞實(shí)驗(yàn)進(jìn)行了比較，得到如下結(jié)論:

(1)利用Fluent動(dòng)網(wǎng)格實(shí)現(xiàn)翼型結(jié)冰后的網(wǎng)格重構(gòu)，可以避免對(duì)單個(gè)結(jié)冰時(shí)間結(jié)束后翼型邊界發(fā)生變化后的流場(chǎng)計(jì)算區(qū)域重新生成網(wǎng)格，提高效率，具有實(shí)際應(yīng)用價(jià)值。

(2)環(huán)境溫度是影響冰形的重要參數(shù)。環(huán)境溫度較高時(shí)，駐點(diǎn)區(qū)域凍結(jié)質(zhì)量小，液態(tài)水向下游溢流，形成具有凹槽的槽狀冰;環(huán)境溫度很低時(shí)，液態(tài)水撞擊到壁面立刻凍結(jié)，外形比較規(guī)則，形成楔形冰;環(huán)境溫度較低時(shí)，冰形兼具槽狀冰和楔形冰的特點(diǎn)，形成混合冰。

(3)利用Messinger結(jié)冰熱力學(xué)模型結(jié)合Fluent動(dòng)網(wǎng)格模塊可以得到和實(shí)驗(yàn)結(jié)果比較吻合的計(jì)算結(jié)果，只是在局部區(qū)域有差異，一定程度上證明了本文所述方法的正確性和可行性，但對(duì)于該模型還需進(jìn)一步改進(jìn)和完善，特別是對(duì)結(jié)冰過(guò)程的機(jī)理需要深入研究。

［1］HAROLD E A，MARK G P，DAVID W S.Modern airfoil ice accretions［R］.AIAA 97-0174，1997.

［2］FORTIN G，PERRON J.CIRAAMIL ice accretion code improvement［R］.AIAA 2009-3968，2009

［3］張大林，陳維建.飛機(jī)機(jī)翼表面霜狀冰結(jié)冰過(guò)程的數(shù)值模擬［J］.航空動(dòng)力學(xué)報(bào)，2003，18(1):88～89.

［4］蔣勝矩，李鳳蔚.基于 N-S方程的翼型結(jié)冰數(shù)值模擬［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2004，22(5):559 ～560.

［5］易賢.飛機(jī)積冰的數(shù)值計(jì)算與積冰試驗(yàn)相似準(zhǔn)則研究［D］.中國(guó)空氣動(dòng)力研究與發(fā)展中心研究生部，2007.

［6］常士楠，艾素霄，等.一種飛機(jī)機(jī)翼表面結(jié)冰過(guò)程仿真方法［J］.系統(tǒng)仿真學(xué)報(bào)，2008，20(10):2538 ～2539.

［7］裘爕綱，韓風(fēng)華.飛機(jī)防冰系統(tǒng)［M］.北京:航空專(zhuān)業(yè)教材審編組，1985.

［8］JEROEN E D，HARRY W M.Numerical simulation of airfoil ice accretion and thermal anti-icing systems［R］.ICAS 2004-7.5.2，2004.

［9］KAYS W M，CRAWFORD M E.Convective heat and mass transfer［M］.Third Edition.Mc Graw-Hill，1993.

［10］JAIWON S，BRIAN B.Prediction of ice shapes and their effect on airfoil drag［J］.Journal of Aircraft，1994，31(2).

［11］DAVID N A，et al.Measurement and correlation of ice accretion roughness［R］.AIAA 98-0486，1998.

［12］VON DOENHOFF，HORTON E A.Low-speed experimental investigation of the effect of sandpaper type of roughness on boundary-layer transition［R］.NACA TN 3858，1956.

［13］MATHIEU F，SAEED F.PARASCHIVOIU I.Surface heat transfer study for ice accretion and anti-icing prediction in three dimension［R］.AIAA 2004-0063，2004.

［14］RUFF G A，BERKOWITZ B M.Users manual for the NASA Lewis ice accretion prediction code(LEWICE)［R］.NASA CR-18519，1990.

［15］MARK G P.LEWICE/E:an Euler based ice accretion code［R］.AIAA 92-0037，1992.

［16］SHIN J，BOND T H.Results of an icing test on a NACA0012 airfoil in the NASA Lewis icing research tunnel［R］.AIAA 92-0647，1992.

［17］MINGIONE G，BRANDI V，ESPOSITO B.Ice accretion prediction on multi-element airfoils［R］.AIAA 1997-0177，1997.