無人機螺旋槳三維積冰計算與積冰控制檢測

栗 樞，曹廣州

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

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

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

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

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

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

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

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

1.1 基本假設

對于旋轉部件，將水膜受到離心力的作用而產生脫離水膜主體的現象稱為“甩脫”。針對“甩脫”現象，做出以下假設。

① 旋轉部件表面的薄水膜流動可以看成不可壓縮的層流流動，且可以假設在小的時間微元內水膜的流動為穩態流動。

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

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

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

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

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

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

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

圖1 積冰控制體

建立積冰與水膜流動相耦合的數學模型，具體如下：

(1)

(2)

(3)

(4)

(5)

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

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

1.3 局部水收集系數的選取

局部水收集系數β是積冰模擬過程中用于表征水滴撞擊特性的無量綱參數，反映出結冰的范圍與程度。局部水收集系數是指當地微元面上的實際水滴撞擊量與相同面積上按自由來流條件獲得的水滴撞擊量之比。采用歐拉-歐拉法計算兩相流場時，可用式(6)和式(7)來計算局部水收集系數。

(6)

(7)

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

1.4 積冰表面水膜甩脫的依據

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

We>2

(8)

得到臨界水膜厚度為

(9)

式中，We為韋伯數，是衡量冰層表面水膜動力學特性的無量綱參數；ρa為氣流密度；Va為當地氣流速度；σiv為表面張力系數。

因此，對螺旋槳水膜甩脫的判據均是基于韋伯數We來判斷的，認為當韋伯數We>2時水膜就會發生甩脫。

1.5 計算方法

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

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

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

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

④ 時間步推進，繼續下一時間步的求解。

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

2.1 模型介紹

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

圖2 扭轉性極強的2714型螺旋槳計算模型

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

圖3 2714型螺旋槳的計算域及網格

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

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

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

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

2.2 網格的劃分和網格獨立性驗證

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

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

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

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

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

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

圖7 不同環境溫度下的冰層厚度分布

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

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可見，隨著水滴含量的增加，結冰位置逐漸向下游發展，這是因為水滴含量增加了，更多的水膜流向下游，邊流動邊凍結成冰，即所謂的明冰。所以在上游的冰型幾乎不變，而在下游處，水滴含量較大的條件使冰型凸起，積冰量變大，逐漸形成更大的角狀冰。

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

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

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

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

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

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

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

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

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

3 結論

本文將適用于旋轉部件表面積冰模擬的模型應用到扭轉性極強的2714型螺旋槳中，介紹了基于2714型螺旋槳強扭轉特性下的模型，并進行了積冰的檢測，通過研究環境溫度分別為-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對螺旋槳表面積冰規律的影響，分析得出以下結論：

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

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

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