常士楠，王 超，趙媛媛，李 延

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

結(jié)冰表面粗糙度預(yù)測與分析

常士楠，王 超，趙媛媛，李 延

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

通過對結(jié)冰表面液態(tài)水的受力分析，確定了水膜、水珠和細流特征粗糙度的大小及其之間的界限劃分，發(fā)展和完善了結(jié)冰表面的粗糙度模型，開發(fā)了耦合該粗糙度模型的多步結(jié)冰計算程序；利用該程序計算結(jié)冰表面的粗糙度大小和分布特征，并對粗糙度影響下的流動和換熱特性進行了計算分析；研究發(fā)現(xiàn)，結(jié)冰表面的粗糙度在分布是不均勻的，且隨著結(jié)冰溫度的降低，局部粗糙度的最大值降低，但在整個結(jié)冰過程中，局部粗糙度的最大值表現(xiàn)出增加的趨勢；粗糙度的存在增強了結(jié)冰表面的對流換熱。

粗糙度；水膜；水珠；細流；結(jié)冰

0 引 言

當(dāng)飛機遭遇含有過冷水滴的結(jié)冰云層時，會發(fā)生結(jié)冰現(xiàn)象。飛機結(jié)冰過程是一個復(fù)雜的氣-液-固三相流動換熱問題，結(jié)冰表面的表面特性如粗糙度等對飛機結(jié)冰表面氣流流動、傳熱與傳質(zhì)、表面撞擊水形態(tài)等產(chǎn)生影響，并最終對結(jié)冰的生長方式和結(jié)冰類型產(chǎn)生重要影響，從而影響因結(jié)冰而導(dǎo)致的氣動性能惡化甚至飛行安全。目前關(guān)于飛機結(jié)冰的研究大多主要基于Messinger在1953年提出的結(jié)冰傳熱模型［1］。該模型揭示了結(jié)冰表面的各項主要傳熱熱流。在給定流場分布和水滴撞擊特性的情況下，利用該模型所建立的能量平衡方程，再結(jié)合結(jié)冰過程中伴隨的質(zhì)量平衡方程，就能夠進行飛機表面的結(jié)冰模擬和預(yù)測。然而這種模擬方法幾乎沒有考慮結(jié)冰表面的物理特性對結(jié)冰的影響，獲得的結(jié)果往往與相同條件下的試驗結(jié)果吻合度不夠理想，很難反映真實的飛機結(jié)冰情況。九十年代初期，Shin等針對結(jié)冰表面的特性發(fā)展了等效粗糙度模型（Shin model）［2］，該模型將結(jié)冰表面粗糙度高度ks表達成結(jié)冰條件的簡單函數(shù)，其計算表達式為：

其中LWC為液態(tài)水含量，g／m3；T∞為環(huán)境溫度，K；u∞為飛行速度，m／s；c為翼型的弦長，m。國內(nèi)外的研究人員大多采用該粗糙度模型開展結(jié)冰模擬。易賢將等效粗糙度模型耦合到結(jié)冰程序中開展了二維和三維結(jié)冰預(yù)測［3］；楊秋明等將等效粗糙度模型和動網(wǎng)格原理結(jié)合，對不同條件下的結(jié)冰展開多步結(jié)冰模擬，并取得了和實驗結(jié)果一定程度的近似［4］。陳維建等在流場計算的湍流模型壁面函數(shù)中加入了粗糙度的影響［5］，但文中關(guān)于粗糙度對流動的影響并非針對飛機結(jié)冰的特殊情況（飛機結(jié)冰表面的雷諾數(shù)遠大于普通流動）。需要指出的是，在給定的結(jié)冰計算條件下等效粗糙度模型的粗糙度高度為一定值，亦即假定機翼整個結(jié)冰表面都具有同樣的粗糙高度，這種假定顯然與實際情況有所區(qū)別。

20世紀80-90年代，Olsen和Hansman利用微攝影和圖像處理技術(shù)對結(jié)冰表面的粗糙度形態(tài)做了大量的研究［6-7］，發(fā)現(xiàn)了結(jié)冰表面三種粗糙度的存在形式：水膜、水珠和細流，并分析了其形成機理；Bragg等研究了粗糙度表面上的邊界層發(fā)展，研究結(jié)果表明大的粗糙度高度導(dǎo)致了翼型表面邊界層由層流到紊流的轉(zhuǎn)變過程比在光滑或小粗糙度翼型表面大大提前［8］；Fortin等建立了粗糙度分析模型，并利用該模型進行了二維結(jié)冰預(yù)測［9］。

盡管基于表面特性的飛機結(jié)冰機理研究已逐漸引起研究人員的重視，但由于影響結(jié)冰的因素非常之多，相互之間又有錯綜復(fù)雜的非線性關(guān)系等，仍然存在許多難點有待突破，如結(jié)冰表面特性（如粗糙度尺度、位置和分布等）的形成規(guī)律、與結(jié)冰條件的關(guān)聯(lián)性等研究少有報道。

本文在總結(jié)以往研究的基礎(chǔ)上，對結(jié)冰表面的粗糙度的形成特點及相關(guān)參數(shù)進行了分析確定，并依據(jù)結(jié)冰表面液態(tài)水的受力特性對表征粗糙度特征的水膜、水珠和細流之間的分界進行了推導(dǎo)界定，進一步發(fā)展了結(jié)冰表面的局部粗糙度模型，并在該粗糙度模型的基礎(chǔ)上研究了結(jié)冰過程中局部粗糙度的大小和分布以及粗糙度影響下的換熱特性，為深入有效的探究飛機結(jié)冰的機理提供理論依據(jù)。

1 局部粗糙度模型的建立

以水膜、水珠和細流為特征的粗糙度大小的確定是通過對壁面上微元控制體進行傳熱、傳質(zhì)分析以及液態(tài)水的受力分析得到的。

1.1 水膜區(qū)粗糙度

水膜表面并非光滑平整，而是存在一定程度的波動，這種波動程度的大小即為水膜區(qū)粗糙度值。取結(jié)冰表面一水膜控制體為研究對象，如圖1所示，其受到的氣動剪切力τi表示為：

式中：μw為水膜的動力粘度，Pa·s；u為速度，m／s；y為液膜的厚度變量，m。式（1）在初始條件y＝0，u＝0下，可得到液膜速度的關(guān)系式為：

則水膜中的平均流速ˉu和水膜的質(zhì)量mw可分別為：

式中ρw為水膜的密度，g／m3；Δb和ef，o分別為液膜控制體的寬度和水膜厚度，m。如圖1所示，對于二維問題，可取Δb＝1。

具體到結(jié)冰過程而言，控制體內(nèi)的液態(tài)水的質(zhì)量

圖1 水膜控制體示意圖Fig.1 Sketch of the film control volume

可由下式求得：

式中f為凍結(jié)系數(shù)，mcap和mevap分別為撞擊到控制體內(nèi)的液態(tài)水量和蒸發(fā)量，kg／s；min為從上一個控制體溢流進的液態(tài)水質(zhì)量，kg／s。

容易得到水膜的厚度為：

根據(jù)水膜區(qū)粗糙度大小的實驗結(jié)果，引入一個校準因子f對水膜厚度進行修正［9］，修正后的水膜厚度為：

式中f＝15，液膜表面的粗糙度可以表示為：

g為重力加速度，m／s2。

由上述各表達式的參數(shù)可看出，水膜區(qū)粗糙度的大小主要取決于水膜受到的剪切力和本身的特性，如水的粘性、密度等，其中對于剪切力τi將在1．3節(jié)給出具體求解過程。

1.2 水珠和細流區(qū)粗糙度

水珠在結(jié)冰表面受到的作用力有氣動力Fd、表面張力Fσ、壁面支持力Fz和重力Fg，如圖2所示。沿氣流方向，水珠在重力切向分量、氣動力、表面張力三者作用下保持穩(wěn)定狀態(tài)，當(dāng)上述力的平衡遭到破壞時，水珠會在氣流的吹襲作用下向后移動，水珠在移動前能夠達到的最大高度eb視為水珠狀粗糙度，其值大小可根據(jù)水珠的受力平衡進行求解。

圖2 水珠受力示意圖Fig.2 Mechanical analysis of the bead on the icing surface

忽略水珠所受到的氣動升力，即認為水珠不會從結(jié)冰表面脫落。此時，水珠受到的氣動力可表示為：

式中ρa為空氣的密度，kg／m3；ue為邊界層外邊界速度，m／s；把水珠等效成一球形顆粒，阻力系數(shù)Cd取值為0．44；Ab為水珠的迎風(fēng)面積，m2，求解式為：

式中：eb為水珠高度，m；θc為水珠未變形時的接觸角度，rad。

水珠的表面張力可以由下面的積分關(guān)系式求得：

σ為水珠的表面張力，N／m；r為水珠底面圓半徑，m；表示為：

式中θ和φ分別為水珠與結(jié)冰表面的接觸角和表面傾斜角度，rad，如圖2所示。在外力作用下水珠與結(jié)冰表面的接觸角θ不斷變化，其變化關(guān)系可認為服從簡單的余弦分布［9］，可表示為：

式中Δθc表示水珠與壁面接觸的滯后角度，rad。

以上各式中涉及到的接觸角θc和滯后角Δθc以及表面張力系數(shù)σ三個參數(shù)的取值大小會直接影響水珠區(qū)粗糙度的計算精度。本文依據(jù)文獻［7］和文獻［10］的試驗結(jié)果，如圖3所示，采用數(shù)值擬合的方法確定各個參數(shù)的表達式如下：

圖3 相關(guān)參數(shù)的擬合Fig.3 Fitting the relevant parameters

上式中te為邊界層溫度，℃。

重力沿結(jié)冰表面切線方向的分量為：

式中V為水珠體積，求解表達式為：

綜合式（10）、式（12）和式（18）可得水珠受力平衡關(guān)系式為：

從以上對于水珠所受的各項力的分析可以看出，獲得結(jié)冰表面邊界層的速度和溫度分布后，即可以通過求解式（20）得到水珠高度eb，對于邊界層的求解見本文第2部分。

假定水珠區(qū)粗糙度值和細流區(qū)粗糙度值相同，大小都為水珠的高度，即：

1.3 粗糙度區(qū)的分界

根據(jù)實驗觀測的結(jié)果，水膜位于結(jié)冰表面的最前緣，水珠位于水膜和撞擊極限之間，細流位于撞擊區(qū)域之外［7］。當(dāng)結(jié)冰條件改變時，三種形態(tài)的分布范圍也會隨之變化。確定不同形態(tài)粗糙度之間的分界對于準確預(yù)測結(jié)冰表面的粗糙度大小具有重要意義。本文通過液膜表面的受力分析對其分界進行確定。

a）水膜與水珠的分界

水膜在沿結(jié)冰表面流動的過程中，受到空氣的剪切應(yīng)力作用，同時，水膜與結(jié)冰表面之間存在著由表面張力引起的內(nèi)聚力，當(dāng)剪切應(yīng)力大于內(nèi)聚力時，水膜能夠保持形態(tài)不變而繼續(xù)向前流動，而當(dāng)剪切應(yīng)力與內(nèi)聚力平衡或者小于內(nèi)聚力時，水膜則會在表面張力的作用下聚集成水珠。因此，水膜與水珠（細流）的分界可通過對結(jié)冰表面液態(tài)水進行受力分析得到。

水膜所受的剪切應(yīng)力的求解表達式為［11］：

式中，u∞為遠場來流速度，m／s；Ai為水膜表面積，m2；Cf為摩擦系數(shù)，其數(shù)值可以通過經(jīng)驗公式得到［12］：

ks表示距離前緣弧長為s處的局部粗糙度，m；稱ks／s為無量綱粗糙度。

由表面張力引起的水膜內(nèi)聚力的計算公式為：

式中，τσ水膜內(nèi)聚力，N；Sσ為水膜周長，m。對于二維問題，微元控制體寬度假定為1，當(dāng)τi≥τσ時，可得：

令Cf，c為臨界摩擦系數(shù)，表示為：

即當(dāng)Cf≥Cf，c時，粗糙度保持水膜狀態(tài)；反之，當(dāng)Cf＜Cf，c時，水膜聚集成水珠。

b）水珠與細流的分界

若水膜在撞擊極限前轉(zhuǎn)變?yōu)樗?，并且在撞擊極限外仍有溢流情況發(fā)生，此時，液態(tài)水會以細流形態(tài)存在。因此，可認為水珠與細流的分界就是水滴撞擊極限。

2 耦合了粗糙度的換熱模型及算法

結(jié)冰翼型壁面的流動非常復(fù)雜，流動存在一個從層流到湍流變化的復(fù)雜過程，轉(zhuǎn)捩點位置的判定對流動類型的區(qū)分十分重要。為此，引入粗糙度雷諾數(shù)的概念［13］。定義

式中，Rek為粗糙度雷諾數(shù)，ν為空氣的運動粘度，m2／s；uk為對應(yīng)粗糙度y＝ks處的空氣速度，m／s，可采用型線假設(shè)法求解［14］：

式中δ為邊界層厚度，m；θlam為層流區(qū)動量邊界層厚度，m。

當(dāng)Rek≤600認為流動為層流，層流區(qū)域中的對流換熱系數(shù)表示成無量綱Stanton數(shù)為：

式中普朗特數(shù)Pr表達式為：

式中Cpa為空氣的比熱容，J／（kg·K）；λ為空氣的導(dǎo)熱系數(shù)，W／（m·K），是邊界層溫度te的函數(shù)：

當(dāng)Rek＞600認為流動為湍流，流動區(qū)域中的對流換熱系數(shù)表示成無量綱Stanton數(shù)為：

式中，Pri為湍流普朗特數(shù)，取值為0.9；Stk為粗糙度Stanton數(shù)，表達式為：

由式（32）和式（33）即可求得湍流區(qū)域的換熱系數(shù)。

由于粗糙度的微觀尺度原因，對粗糙度模型的預(yù)測精度提出了較高的要求。為了有效的對所建立的局部粗糙度模型進行評估，實現(xiàn)對結(jié)冰表面粗糙度大小和分布特點及粗糙度影響下的換熱特性的精確預(yù)測，開發(fā)了耦合局部粗糙度模型的多步結(jié)冰預(yù)測程序，算法流程如圖4所示。

圖4 計算流程圖Fig.4 Flow chart of computation

3 算例分析

選取NACA0012翼型為計算模型，并以結(jié)冰時間為360s時的計算結(jié)果分別從結(jié)冰表面平均粗糙度和冰形兩方面與文獻數(shù)據(jù)進行對比分析，其中，流場計算部分借助Fluent計算流體軟件［15］完成，結(jié)冰模型采用經(jīng)典的Messinger熱質(zhì)耦合模型［1］。計算條件列于表1中。

表1 計算條件Table 1 Calculation conditions

3.1 平均粗糙度預(yù)測

針對局部粗糙分析模型，引入平均沙粒粗糙度的概念，即對結(jié)冰表面上所有微元控制體內(nèi)的粗糙度進行積分，求出積分平均值，該值即為等效沙粒粗糙高度，表示為式（34）。選取的計算溫度為-4．4℃，-6．1℃，-7．8℃，-10℃，-13．35℃，-19．4℃，-28．4℃。

式中ds為微元控制體弧段長度，m。

圖5給出了各溫度條件下，采用局部粗糙度模型所得的平均粗糙度值與Shin等效粗糙度模型和Fortin模型的對比結(jié)果。本文開始部分已經(jīng)對等效粗糙度模型進行了說明，由式（1）也可以看出，在給定結(jié)冰條件下，由該模型計算得到的粗糙高度與溫度呈線性關(guān)系。文中的預(yù)測結(jié)果與Fortin模型的預(yù)測值較為接近。局部粗糙度模型是基于壁面微元控制體內(nèi)的熱質(zhì)平衡和對液態(tài)水的受力分析建立起來的模型，綜合考慮了環(huán)境溫度、氣流速度、液態(tài)水含量（LWC）等結(jié)冰條件以及液態(tài)水與固體壁面的接觸角度、長度和面積等多種幾何因素的影響，因此，導(dǎo)致了采用局部粗糙度模型得到平均粗糙度值具有較高程度的非線性分布。

圖5 不同溫度時平均粗糙度大小對比Fig.5 Comparison of average roughness at different temperature

3.2 冰形預(yù)測

圖6 混合冰結(jié)冰條件下的冰形比較Fig.6 Comparison of ice shapes under the mixedicing condition

冰形是耦合了局部粗糙度模型的多步結(jié)冰程序的最終計算結(jié)果，是考量模型算法的重要標(biāo)準。從圖6中可以看出，本文的計算程序?qū)τ谛螤顝?fù)雜的混合冰的冰角位置、結(jié)冰冰形以及結(jié)冰量等預(yù)測結(jié)果與實驗結(jié)果吻合較好，表明文中耦合了局部粗糙度模型的算法程序用于過冷水滴的撞擊結(jié)冰是成功的。

4 耦合局部粗糙度模型的計算結(jié)果與分析

通過對局部粗糙度模型的分析可知，水膜狀、水珠狀以及細流狀粗糙度的分布特征是結(jié)冰氣象條件和流體力學(xué)因素共同作用的結(jié)果，因此，探究不同條件下的結(jié)冰表面的粗糙度分布范圍和大小，對于從微觀角度揭示飛機結(jié)冰機理具有重要意義。

4.1 結(jié)冰過程中局部粗糙度大小與分布

圖7 不同溫度條件下結(jié)冰表面的局部粗糙度分布Fig.7 Local roughness distribution on the icing surface for different temperatures

圖7給出了結(jié)冰過程中結(jié)冰表面局部粗糙度大小和分布特點。結(jié)冰過程中，由于結(jié)冰表面液態(tài)水形態(tài)的變化，促使冰層的不斷累積，從而造成結(jié)冰表面流場的不斷變化，而流場的變化又會影響結(jié)冰表面液態(tài)水的量與分布，從而引起局部粗糙度的大小和分布變化。容易看出，利用本文程序預(yù)測的粗糙度在結(jié)冰表面上呈“M”形分布，駐點處的粗糙度值在整個結(jié)冰過程中都是最小的；沿弦長氣流運動方向，粗糙度值逐漸增大，并在一定位置處達到最大值；在結(jié)冰過程中，粗糙度的最大值表現(xiàn)出增加的趨勢，但是綜合三種溫度條件下的粗糙度分布而言，粗糙度的最大值（-4．4℃，ks，max＝2．72mm；-6．1℃，ks，max＝2．45mm；-10℃，ks，max＝1．6mm）是隨著環(huán)境溫度的降低而降低的，這是因為溫度的降低促使結(jié)冰表面的液態(tài)水量的減少，從而使局部粗糙度變小。

4.2 結(jié)冰表面對流換熱特性分析

圖8 結(jié)冰過程中結(jié)冰表面的Rek分布Fig.8 Distribution of Rekon the surface during the icing process

結(jié)冰表面氣流的流動形態(tài)對表面的對流換熱起著重要作用，其中層流區(qū)與湍流區(qū)的分布和范圍大小反應(yīng)了結(jié)冰表面的對流換熱特性。圖8給出了不同環(huán)境溫度條件下結(jié)冰表面的粗糙度雷諾數(shù)分布的計算結(jié)果。粗糙度雷諾數(shù)Rek為600的水平線把Rek分布曲線與橫向坐標(biāo)軸圍成的區(qū)域分成了上下兩部分，水平線上方區(qū)段對應(yīng)的弧長范圍表示結(jié)冰表面的湍流區(qū)域，水平線下方區(qū)段對應(yīng)的弧長范圍表示結(jié)冰表面的層流區(qū)域。從圖中可以看出，特定溫度下的結(jié)冰過程中，結(jié)冰表面的流動形態(tài)的區(qū)域和大小幾乎沒有發(fā)生變化；隨著溫度的降低（-4．4℃，-6．1℃，-10℃），結(jié)冰表面的湍流區(qū)域范圍減小，且表征流動形態(tài)的粗糙度雷諾數(shù)的最大值也減小。對比圖7可以看出，湍流區(qū)對應(yīng)著粗糙度值較大的分布區(qū)域。

結(jié)冰表面流動區(qū)域內(nèi)的對流換熱系數(shù)大小可用無量綱對流換熱系數(shù)Stanton數(shù)（St）來表示。當(dāng)流動為層流時，St與粗糙度的關(guān)系如圖9（a）所示，St隨著無量綱粗糙度ks／s的增加而增加，因此，當(dāng)同一弧長位置處的粗糙度增加時，St增加，即結(jié)冰表面的對流換熱作用增強；當(dāng)流動為湍流形態(tài)時，St主要受粗糙度、粗糙度Stanton數(shù)（Stk）的影響，其影響關(guān)系如圖9（b）所示。無量綱粗糙度ks／s對St的影響作用與流動為層流時相同；當(dāng)ks／s不變時，St隨著Stk的增加表現(xiàn)出增加的趨勢，但是當(dāng)Stk＞10時，St不再隨Stk的增加而發(fā)生變化，即此時可以忽略粗糙度Stanton數(shù)Stk對結(jié)冰表面的對流換熱作用的影響。

從以上的分析可知，相同的條件下結(jié)冰表面的粗糙度越大，其表面的對流換熱作用越強。

圖9 無量綱粗糙度ks／s對St的影響Fig.9 Effect of ks／s on St

5 結(jié) 論

本文在對結(jié)冰表面液態(tài)水的受力分析的基礎(chǔ)上發(fā)展了局部粗糙度模型，開發(fā)了耦合局部粗糙度模型的多步結(jié)冰程序，并對結(jié)冰表面的粗糙度分布特性和粗糙度影響下的結(jié)冰表面的傳熱特性進行了研究，得出以下主要結(jié)論：

（1）局部粗糙度模型在一定程度上反映了結(jié)冰表面粗糙度的分布特點。利用耦合了該粗糙度模型的多步結(jié)冰程序?qū)Y(jié)冰表面粗糙度和結(jié)冰冰形進行了模擬計算，計算結(jié)果有助于從微觀角度認識過冷水滴的撞擊結(jié)冰機理。

（2）結(jié)冰表面的局部粗糙度的分布是不均勻的，但在整個結(jié)冰過程中粗糙度的整體分布是穩(wěn)定的；在本文的計算條件下，粗糙度較大的區(qū)域氣流流動的形態(tài)由層流轉(zhuǎn)變?yōu)橥牧?，粗糙度促進了流動區(qū)域內(nèi)的對流換熱。

（3）局部粗糙度模型綜合考慮了結(jié)冰表面不同區(qū)域粗糙度的大小和形態(tài)，與實際的結(jié)冰表面比較近似，在今后的結(jié)冰代碼的改進與開發(fā)中應(yīng)著重考慮局部粗糙度模型。

［1］MESSINGER B L．Equilibrium temperature of an unheated icing surface as a function of airspeed［J］．Journal of the Aeronautical Sciences，1953，20（1）：29-42．doi：10．2514／8．2520

［2］JAIWON Shin，BRIAN Berkowitz，HSUN Chen，et al．Prediction of ice shapes and their effect on airfoil performance［R］．AIAA 91-0264．

［3］YI X．Numerical computation of aircraft icing and study on icing test scaling law［D］．［PHD Thesis］．China Aerodynamics Research and Development Center，2007．（in Chinese）

易賢．飛機積冰的數(shù)值計算與積冰試驗相似準則研究［D］．［博士學(xué)位論文］．綿陽：中國空氣動力研究與發(fā)展中心，2007．

［4］CHANG S N，YANG Q M，LI Y．Quasi-steady numerical simulation of ice accretion on airfoil［J］．ACTA Aerodynamica Sinica，2011，29（06）：302-308．（in Chinese）

常士楠，楊秋明，李延．翼型表面結(jié)冰準定常數(shù)值模擬［J］．空氣動力學(xué)學(xué)報，2011，29（06）：302-308．

［5］CHEN W J，ZHANG D L．Numerical simulation of glaze ice accretion process［J］．Journal of Aerospace Power，2005，20（3）：472-476．（in Chinese）

陳維建，張大林．瘤狀冰結(jié)冰過程的數(shù)值模擬［J］．航空動力學(xué)報，2005，20（3）：472-476．

［6］OLSEN W，WALKER E．Experimental evidence for modifying the current physical model for ice accretion on aircraft surfaces［R］．NASA TM-87184，1986．

［7］HANSMAN R J Jr，REEHORST A，SIMS J．Analysis of surface roughness generation in aircraft ice accretion［R］．AIAA 92-0298．

［8］BRAGG M B，CUMMINGS M J，HENZE C M．Boundary-layer and heat-transfer measurements on an airfoil with simulated ice roughness［R］．AIAA 96-0866．

［9］GUY F．Prediction of 2D airfoil ice accretion by bisection method and by rivulets and beads modeling［R］．AIAA 2003-1076．

［10］LUXFORD G．Experimental and modelling investigation of the deformation，drag and break-up of drizzle droplets subjected to strong aerodynamics forces in relation to SLD aircraft icing［D］．Cranfield university school of engineering，2005．

［11］TOMPSON B E，MARROCHELLO M R．Rivulet formation in surface-water flow on an airfoil in rain［J］．AIAA Journal，1999，37（1）：45-49．

［12］BEAUGENDRE H，MORENCY F，HABASHI W G．Roughness implementation in fensap-ice：model calibration and influence on ice shapes［J］．Journal of Aircraft，2003，40（6）：1212-1215．doi：10．2514／2．7214

［13］Von DONEHOFF A，HORTON E．A low-speed experimental investigation of the effect of a sandpaper type of roughness on boundary-layer transition［R］．NACA Report No．1349，1956．

［14］ZHANG Z X，DONG Z N．Viscous fluid mechanics［M］．Beijing：Tsinghua University Press，1999．（in Chinese）

章梓雄，董曾南．粘性流體力學(xué)［M］．北京：清華大學(xué)出版社，1999．

［15］FLUENT 6．3 User′s Guide［M］．2006．

Modeling of roughness dimension and distribution on the icing surface

CHANG Shinan，WANG Chao，ZHAO Yuanyuan，LI Yan
（School of Aeronautic Science and Engineering，Beihang University，Beijing 100191，China）

It′s recognized that the roughness dimension and distribution on the icing surface influence the flow structure in the boundary layer and thus the heat transfer characteristics is affected further．Therefore，constructing an appropriate roughness model becomes a strong requirement in order to predict the ice shape accurately．In the current work，a local roughness model was developed based on the mechanical analysis of the liquid water on the icing surface．The proposed roughness model is characterized by film，beads and rivulet．The boundaries between film，beads and rivulet are defined according to the experimental observations in references．The set of transient governing equations developed were solved numerically with the finite vol-ume method using a computational fluid dynamics（CFD）code．The mathematical model can be used to analyze the effect of different icing conditions on the roughness and thus the characteristics of the heat transfer in the icing process can be obtained．Numerical results are presented to show the applicability of the developed model and can be helpful to understand the icing mechanism further．In addition，the multi-step ice predicting code developed in this work can be applied in anti／deicing system design directly．The predictions show that the roughness distribution on icing surface is non-uniform，and the maximum roughness is decreasing with temperature decreases．However，during the whole icing process，the maximum roughness performs an increasing trend．The convective heat transfer was enhanced due to the roughness on icing surface．

roughness；film；beads；rivulet；ice accretion

O359

Adoi:10.7638／kqdlxxb-2012.0177

0258-1825（2014）05-0660-08

2012-10-30；

2013-01-01

國家自然科學(xué)基金（10872020，11072019，11372026）

常士楠（1968-），教授，研究方向：飛機和發(fā)動機結(jié)冰及防／除冰研究．E-mail：sn＿chang＠buaa．edu．cn

常士楠，王 超，趙媛媛，等．結(jié)冰表面粗糙度預(yù)測與分析［J］．空氣動力學(xué)學(xué)報，2014，32（5）：660-667．

10．7638／kqdlxxb-2012．0177． CHANG S N，WANG C，ZHAO Y Y，et al．Modeling of roughness dimension and distribution on the icing surface［J］．ACTA Aerodynamica Sinica，2014，32（5）：660-667．