栗 樞，曹廣州

(1.南京航空航天大學(xué) 能源與動力學(xué)院，江蘇 南京 210016；2.南京航空航天大學(xué) 無人機研究院，江蘇 南京 210016)

當(dāng)飛機穿過含有過冷水滴的云層時，空氣中的過冷水滴會撞到發(fā)動機進口、機翼無人機螺旋槳等迎風(fēng)部件表面，并產(chǎn)生結(jié)冰現(xiàn)象。這些迎風(fēng)部件表面形成的結(jié)冰現(xiàn)象會破壞表面氣動外形，危及飛機安全，甚至發(fā)生墜毀。因此研究飛機積冰問題，在工程上具有十分重要的意義。

能量有限、結(jié)構(gòu)受限、對重量敏感等特點導(dǎo)致中小型無人機的防除冰設(shè)計相對困難，效費比太低，因此，在有人機中成熟使用的防除冰系統(tǒng)一般不能直接用于無人機上。目前來看，用戶們普遍采取的應(yīng)對方法是在環(huán)境條件較好的時間段進行無人機的飛行，比如：科索沃戰(zhàn)爭中美軍為避免結(jié)冰危害，不得不在10月至4月期間禁飛無人機。這樣就會大大影響無人機的出勤率。即便如此，仍有25%的架次遭遇結(jié)冰，導(dǎo)致不同程度的故障或墜毀。所以說，結(jié)冰是不可避免的，即使是環(huán)境條件相對好的情況下會產(chǎn)生結(jié)冰現(xiàn)象。

國外對有人機迎風(fēng)表面積冰的數(shù)值模擬在20世紀40年代就已開展，1953年Messinger首次提出用于積冰數(shù)值計算的模型，此后該模型得到了廣泛的應(yīng)用。20世紀70年代，隨著CFD的發(fā)展，積冰數(shù)值模擬技術(shù)快速發(fā)展，80年代之后陸續(xù)出現(xiàn)了一批基于Messinger模型開發(fā)的積冰模擬/防冰設(shè)計軟件，如美國的LEWICE[1]、英國的 TRAJICE[2]、法國的 ONERA[3]、意大利的 CIRA[4]以及加拿大的 FENSAP-ICE[5]、Myers[6-7]等。但Messinger模型只是對霜冰的計算分析較為準確，認為當(dāng)前控制體中未凝結(jié)的液態(tài)水全部流入下一個控制體中，沒有考慮冰層表面薄水膜流動對積冰的影響[2]。

國內(nèi)絕大多數(shù)的積冰模型都是以Messinger模型為基礎(chǔ)開發(fā)的二維積冰模型，而曹廣州[8-9]則在處理水膜流動的問題上通過量綱分析的方法簡化了N-S方程，得到了3個水膜流動的方程，建立了國內(nèi)第一個考慮水膜流動的三維積冰模型并開發(fā)其計算方法[10]。該模型認為水膜在重力和空氣剪切力的共同作用下會分別向展向和弦向流動，并對國外典型的三維積冰進行了數(shù)值模擬研究，結(jié)果匹配良好。

對于旋轉(zhuǎn)部件，陳寧立[11-12]在曹廣州普通迎風(fēng)部件表面三維積冰與薄水膜流動耦合的基礎(chǔ)上，考慮了離心力對水膜流動的影響，發(fā)展了旋轉(zhuǎn)部件表面薄水膜流動與三維積冰相變耦合的數(shù)學(xué)模型，并對直升機螺旋槳、發(fā)動機進口整流罩等典型旋翼結(jié)構(gòu)進行了數(shù)值模擬研究，又通過與實驗結(jié)果對比驗證了考慮離心力作用下該三維積冰模型的準確性。梁鵬[13-14]自主搭建小冰風(fēng)洞，通過可視化方法得到不同條件下的積冰形狀，并通過圖像分析了不同條件下的結(jié)冰規(guī)律。

本文考慮到螺旋槳翼型的三維扭轉(zhuǎn)特性，開發(fā)了適合螺旋槳三維積冰的結(jié)構(gòu)化-非結(jié)構(gòu)化的混合網(wǎng)格劃分方法。通過CFX流場計算，找到幾組最佳速度-轉(zhuǎn)速匹配，在考慮到飛行來流和旋轉(zhuǎn)來流情況下使三維積冰模型更貼合螺旋槳氣動外形。又考慮到旋轉(zhuǎn)條件下水膜的甩脫，提出考慮飛行來流的旋轉(zhuǎn)條件下的撞擊特性和三維積冰模型的數(shù)學(xué)模型及其計算方法，利用Fortran開發(fā)了相應(yīng)的積冰模擬程序，對整個積冰過程進行檢測，并進行了不同速度/轉(zhuǎn)速、溫度、水滴含量和水滴直徑下的積冰數(shù)值模擬和對比分析。

1 螺旋槳表面三維積冰模型與檢測方法介紹

1.1 基本假設(shè)

對于旋轉(zhuǎn)部件，將水膜受到離心力的作用而產(chǎn)生脫離水膜主體的現(xiàn)象稱為“甩脫”。針對“甩脫”現(xiàn)象，做出以下假設(shè)。

① 旋轉(zhuǎn)部件表面的薄水膜流動可以看成不可壓縮的層流流動，且可以假設(shè)在小的時間微元內(nèi)水膜的流動為穩(wěn)態(tài)流動。

② 普通積冰部件表面未凝結(jié)的水會形成連續(xù)的薄水膜，該水膜厚度的量級比較小[15]，一般在10-4m量級，旋轉(zhuǎn)部件考慮到離心力的作用，并考慮到轉(zhuǎn)速的不固定性，水膜厚度的量級一般在10-7，轉(zhuǎn)速相對較高的無人機螺旋槳可達10-8。

③ 暫認為水膜主體的能量、動量的不受水膜甩脫的影響。

④ 暫認為周圍的氣液稀疏兩相流場不受水膜甩脫的影響。

⑤ 暫不考慮水滴撞擊所帶來的濺灑。

⑥ 旋轉(zhuǎn)積冰部件表面的薄水膜的溫度一般很低，因此暫不考慮輻射換熱所帶來的影響。

1.2 螺旋槳表面水膜流動與三維積冰模型

積冰的形成首先是水滴撞擊到飛機的迎風(fēng)部件表面形成水膜，之后在空氣-過冷水滴兩相流場中對流換熱逐漸形成冰層。所采用的積冰控制體如圖1所示，該控制體由3個部分組成，最上方的是空氣-過冷水滴兩相流場，中間是未凍結(jié)水膜，最下方則是冰層。其中，Hw、Tw分別為水膜厚度和溫度；Hi、Ti分別為冰層厚度和溫度；mimp為水滴撞擊的質(zhì)量；mice為積冰的質(zhì)量；min為流入該控制體的質(zhì)量；mout為流向下一個控制體的質(zhì)量。

圖1 積冰控制體

建立積冰與水膜流動相耦合的數(shù)學(xué)模型，具體如下：

(1)

(2)

(3)

(4)

(5)

式中，ρi、ρw分別為冰和水的密度；rx、ry、rz分別為r在x、y、z方向上的分量，r為該點到旋轉(zhuǎn)軸的距離向量的反向;Lf為積冰的相變潛熱，本文取334400 J/kg；λi、λw分別為冰和水的導(dǎo)熱系數(shù)。

其中，式(1)～式(3)分別為冰層表面水膜流動的連續(xù)方程、動量方程和能量方程，式(4)為冰層中的能量方程，式(5)為水膜的積冰相變方程。

1.3 局部水收集系數(shù)的選取

局部水收集系數(shù)β是積冰模擬過程中用于表征水滴撞擊特性的無量綱參數(shù)，反映出結(jié)冰的范圍與程度。局部水收集系數(shù)是指當(dāng)?shù)匚⒃嫔系膶嶋H水滴撞擊量與相同面積上按自由來流條件獲得的水滴撞擊量之比。采用歐拉-歐拉法計算兩相流場時，可用式(6)和式(7)來計算局部水收集系數(shù)。

(6)

(7)

式中，αw為壁面上水滴的體積含量；Unw為壁面上水滴的法向撞擊速度，取壁面外法線為正方向，因此Unw均為負值，式(6)中的負號是為了保證局部水收集系數(shù)β的值均為正；U0為螺旋槳表面受到合來流速度，數(shù)值上等于飛行來流U∞和旋轉(zhuǎn)來流Ur的合速度，相比于普通翼型的水收集系數(shù)計算多了一項旋轉(zhuǎn)來流。

1.4 積冰表面水膜甩脫的依據(jù)

當(dāng)氣流曳力做的功和水膜由于受到離心力增加的潛能大于水膜和螺旋槳表面的黏附功時，水膜會發(fā)生甩脫。具體得出：

We>2

(8)

得到臨界水膜厚度為

(9)

式中，We為韋伯?dāng)?shù)，是衡量冰層表面水膜動力學(xué)特性的無量綱參數(shù)；ρa為氣流密度；Va為當(dāng)?shù)貧饬魉俣?；σiv為表面張力系數(shù)。

因此，對螺旋槳水膜甩脫的判據(jù)均是基于韋伯?dāng)?shù)We來判斷的，認為當(dāng)韋伯?dāng)?shù)We>2時水膜就會發(fā)生甩脫。

1.5 計算方法

本文的計算方法步驟如下：

① 積冰從未積冰的光滑壁面開始，t=0時，Hw=0，Hi=0，開始進行三維建模、網(wǎng)格劃分、兩相流場計算，并通過自開發(fā)軟件提取流場數(shù)據(jù)。

② 計算得到積冰厚度、水膜厚度等參數(shù)。

③ 計算We>2是否成立，若滿足，計算求解新的水膜厚度。

④ 時間步推進，繼續(xù)下一時間步的求解。

2 螺旋槳表面三維積冰過程的檢測內(nèi)容與積冰控制

2.1 模型介紹

選取2714型螺旋槳進行三維積冰模擬，如圖2所示，槳的直徑27 in，約為0.6858 m(1 in=0.0254 m)。由于前方機身的遮擋，螺旋槳轉(zhuǎn)軸附近受流場影響較弱，故做簡化處理，去除螺旋槳轉(zhuǎn)軸附近的槳葉部分，由于槳尖效率低且結(jié)構(gòu)復(fù)雜，計算模型半徑縮減至為0.31 m。

圖2 扭轉(zhuǎn)性極強的2714型螺旋槳計算模型

由于中小型無人機一般采用兩葉螺旋槳，故只選取一半進行研究，所建立的計算域為圖3所示的半圓柱結(jié)構(gòu)，圖3(a)中紅色區(qū)域為整個計算域，藍色區(qū)域為模擬的機身部分，黃色區(qū)域為內(nèi)部結(jié)構(gòu)化網(wǎng)格部分，包含螺旋槳，同時考慮前面機身遮擋的影響而挖去一部分。螺旋槳的弦長在各個截面處都不相同，范圍為0.03～0.05 m之間，內(nèi)部計算域為提高計算精度，螺旋槳所受到的旋轉(zhuǎn)來流的上游和下游取單倍弦長，厚度取槳長。而在外部非結(jié)構(gòu)化區(qū)域，要考慮部分機身的長度，飛行來流的上游和下游為2∶1的15倍弦長，圓柱的半徑為1.5倍的槳長。半圓柱的兩個側(cè)面采用周期性邊界條件，圖3(b)為內(nèi)部結(jié)構(gòu)化網(wǎng)格的放大圖。利用ANSYS-ICEM對螺旋槳進行網(wǎng)格劃分，為了減少網(wǎng)格量，同時確保計算精度，采用了結(jié)構(gòu)-非結(jié)構(gòu)的混合網(wǎng)格，如圖3所示，對螺旋槳的近壁區(qū)采用質(zhì)量較好、密度較大的結(jié)構(gòu)化網(wǎng)格。而在外部區(qū)域采用較為粗糙的非結(jié)構(gòu)化網(wǎng)格。采用表1的參數(shù)進行流場設(shè)置和數(shù)據(jù)提取。

圖3 2714型螺旋槳的計算域及網(wǎng)格

表1 2714型螺旋槳的計算條件

圖4為選取不同半徑位置螺旋槳的截面圖，以螺旋槳各個截面上最左端為駐點位置，上表面即壓力面的弧長坐標(biāo)值為正值，下表面即吸力面坐標(biāo)值為負值。由圖4可見，螺旋槳在不同的半徑位置的截面在形狀、角度甚至幾何尺度上都各不相同，體現(xiàn)了其強扭轉(zhuǎn)特性。由1.3節(jié)局部水收集系數(shù)公式可知，其本身與氣流速度有關(guān)，這正與要研究的流場速度-轉(zhuǎn)速匹配相關(guān)聯(lián)，其次本文最終旨在研究積冰的影響，局部水收集系數(shù)作為積冰的重要參數(shù)，影響著整個積冰過程和結(jié)果，故選用該參數(shù)作為一些算例展示及驗證分析。

螺旋槳表面的局部水收集系數(shù)分布曲線如圖4所示，其中橫坐標(biāo)為螺旋槳的弧長坐標(biāo)。由圖4可知，壓力面即S>0方向的水滴撞擊極限要遠遠大于吸力面即S<0方向上的水滴撞擊極限，因此局部水收集系數(shù)呈不對稱分布。且在S>0方向，旋轉(zhuǎn)半徑越大的位置，水滴撞擊極限越廣。這是由于隨著旋轉(zhuǎn)半徑的增大，使得線速度增大，來流水滴撞擊速度增大，水滴慣性也隨之增大，水滴就更容易撞擊到螺旋槳表面，其局部水收集系數(shù)也相應(yīng)增大。

圖4 光滑螺旋槳不同半徑位置截面形狀圖

2.2 網(wǎng)格的劃分和網(wǎng)格獨立性驗證

在算例1的條件下，通過網(wǎng)格單元數(shù)分別為65萬、85萬、100萬、120萬、140萬進行流場計算，并通過數(shù)據(jù)提取得到相應(yīng)的局部水收集系數(shù)數(shù)據(jù)，如圖5所示，工程上螺旋槳總半徑的0.7倍的位置一般為螺旋槳的最初設(shè)計點。故取螺旋槳表面在r/R=0.7截面上的局部水收集系數(shù)，由圖5可知，當(dāng)網(wǎng)格單元數(shù)量達到100萬時，局部水收集系數(shù)基本上不再變化，可以判斷在所采用的網(wǎng)格疏密范圍內(nèi)，計算結(jié)果受網(wǎng)格量的影響可以忽略不計。最終采用網(wǎng)格單元數(shù)為100萬的網(wǎng)格進行計算分析。采取網(wǎng)格單元的弦長隨著展向位置的不同而不斷變化，故本文在不加說明的情況下，均取螺旋槳總長度的0.7倍位置處的積冰參數(shù)進行研究。

圖5 網(wǎng)格獨立性驗證(r/R=0.7)

2.3 不同溫度對螺旋槳積冰影響的檢測

選取溫度為-5 ℃、-10 ℃、-15 ℃、-20 ℃，研究不同溫度對螺旋槳積冰的影響。圖6為T=-10 ℃下r/R=0.7截面處的局部水收集系數(shù)分布。由局部水收集系數(shù)公式可知，局部水收集系數(shù)僅與水滴密度、體積分數(shù)、水滴含量、速度有關(guān)，與環(huán)境溫度的高低無關(guān)，而且溫度的改變不會影響水滴撞擊量的改變，所以4種不同環(huán)境溫度下的局部水收集系數(shù)分布相同。

圖6 T=-10 ℃下局部水收集系數(shù)分布

不同環(huán)境溫度下r/R=0.7截面位置冰層厚度分布情況如圖7所示。由圖7可見，隨著環(huán)境溫度的降低，相同截面上的冰形高度越大，冰形的覆蓋范圍也越大。

圖7 不同環(huán)境溫度下的冰層厚度分布

不同環(huán)境溫度下r/R=0.7截面處的最終冰型輪廓圖如圖8所示。由圖8可知，隨著溫度的降低，結(jié)冰位置基本上不受影響，因為未能從氣動上改變整個流場，只是溫度更低了，水滴更易凝結(jié)成冰，在相同的位置上，積冰量更大，冰型更厚，逐漸形成更大的角狀冰。當(dāng)水滴撞擊量相同時，溫度越低，溫度的變化量增大，水膜表面與兩相流場間的換熱量增大，相同質(zhì)量流量下水滴吸收的熱量增大，根據(jù)熱力學(xué)平衡定律，螺旋槳表面的積冰量就變大。同時環(huán)境溫度越低，空氣能吸收的熱量也越多，由熱力學(xué)平衡定律可知，相應(yīng)的積冰放熱量增大，螺旋槳表面的積冰量增大。

2.4 不同水滴含量對螺旋槳積冰影響的檢測

選取水滴含量分別為0.5 g/m3、1.0 g/m3、1.5 g/m3、2.0 g/m3,研究不同水滴含量對螺旋槳積冰的影響。不同水滴含量下r/R=0.7截面位置冰層厚度和最終冰型輪廓圖分布情況如圖9、圖10所示。由圖9可見，隨著水滴含量的增加，結(jié)冰位置逐漸向下游發(fā)展，這是因為水滴含量增加了，更多的水膜流向下游，邊流動邊凍結(jié)成冰，即所謂的明冰。所以在上游的冰型幾乎不變，而在下游處，水滴含量較大的條件使冰型凸起，積冰量變大，逐漸形成更大的角狀冰。

圖9 不同水滴含量的冰層厚度分布

圖10 不同水滴含量的冰型輪廓圖

2.5 不同水滴直徑對螺旋槳積冰影響的檢測

選取水滴直徑分別為15 μm、20 μm、25 μm、35 μm，研究不同水滴直徑對螺旋槳表面積冰的影響，不同水滴直徑下r/R=0.7截面的局部水收集系數(shù)分布情況如圖11所示。由圖11可知，隨著水滴直徑的增大，局部水收集系數(shù)有所增加。因為水滴直徑越大，水滴的慣性就越大，從而使水滴更容易撞擊到螺旋槳表面。

圖11 不同水滴直徑的局部水收集系數(shù)分布

不同水滴直徑下r/R=0.7截面的冰層厚度分布情況如圖12所示。由圖12可知，隨著水滴直徑的增加，積冰的范圍和厚度均有所擴大和增加。

圖12 不同水滴直徑的冰層厚度分布

不同水滴直徑下r/R=0.7截面的最終冰型輪廓圖如圖13所示。由圖13可知，上游的冰型輪廓基本保持一致，在下游處，隨著水滴直徑的增加，冰型有所擴大。

圖13 不同水滴直徑的冰型輪廓圖

3 結(jié)論

本文將適用于旋轉(zhuǎn)部件表面積冰模擬的模型應(yīng)用到扭轉(zhuǎn)性極強的2714型螺旋槳中，介紹了基于2714型螺旋槳強扭轉(zhuǎn)特性下的模型，并進行了積冰的檢測，通過研究環(huán)境溫度分別為-5 ℃、-10 ℃、-15 ℃、-20 ℃，水滴含量分別為0.5 g/m3、1.0 g/m3、1.5 g/m3、2.0 g/m3以及水滴直徑分別為15 μm、20 μm、25 μm、35 μm對螺旋槳表面積冰規(guī)律的影響，分析得出以下結(jié)論：

① 隨著環(huán)境溫度的降低，局部水收集系數(shù)和對流換熱系數(shù)基本不受影響，檢測出空氣中的水滴更容易結(jié)冰，冰層厚度增大，水膜厚度增大，冰型輪廓逐漸擴大。

② 隨著水滴含量的增大，檢測出冰層厚度增大，水膜厚度增大，上游的冰型輪廓幾乎不變，更多的水滴流向下游，從而形成明顯的角狀冰。

③ 隨著水滴直徑的增大，局部水收集系數(shù)增大，對流換熱系數(shù)幾乎不受影響，檢測出積冰量增大，水膜厚度略有增加，上游積冰形態(tài)基本上不變，下游的積冰量越來越大，逐漸形成明顯的角狀冰。