葉型結冰計算及流場分析

任永鵬，申連洋，王忠義，王艷華，萬雷，王松

(1. 哈爾濱工程大學 動力與能源工程學院， 哈爾濱 150001)

(2. 海裝沈陽局駐哈爾濱地區第二軍事代表室， 哈爾濱 150001)

0 引 言

飛機機翼、壓氣機進氣部件結冰主要是由于過冷水滴(0 ℃以下仍為液態的水滴)撞擊到固體表面冷凝，因此產生了結冰現象[1]。大氣環境溫度、液態水含量、過冷液滴直徑和氣流速度等因素都會對結冰的種類、結冰量和結冰區域產生影響[2-3]。翼型和葉型結冰都會改變其原有的氣動外形，影響氣動特性。因此，對于結冰機理的研究十分重要[4-5]。

20世紀60年代，科研人員開始對結冰問題開展數值模擬研究。21世紀隨著計算機科學的迅猛發展，數值模擬已成為預測翼型結冰的重要手段之一。以計算流體力學和計算傳熱學為基礎，流場、液滴撞擊特性、結冰量的計算方法日趨成熟，數值模擬以其經濟高效的優勢，在結冰問題研究上得到了廣泛的應用。

國外，美國代頓大學的C.MacArthur[6]建立了二維翼型上明冰和霜冰積累的數學模型，對含有液態水的結冰過程進行熱力學分析，使液滴收集系數、傳熱和物質交換隨著翼型表面積冰形狀的改變而進行動態更新，從而對翼型周邊的流場及液滴軌跡進行標準化的計算；美國通用電氣全球研究中心的S.Saxena等[7]，模擬了液滴在動葉、靜葉上的飛濺效應、葉片固結效應和葉片尾緣的蔓延效應，研究發現液滴與動、靜葉間的相互作用會導致壓氣機喘振裕度的變化達到25%；S.Ozgen等[8]利用FORTRAN代碼對二維NACA0012翼型積冰形狀和二維Twin Otter翼型的水滴收集系數進行預測，計算結果與試驗結果吻合較好。我國對于翼型結冰的研究雖起步較晚，但近些年來在結冰數值模擬方法上也積累了一定的科研成果，趙秋月[9]利用Fluent用戶自定義設計，計算過冷水滴的運動軌跡、航空發動機三維進口支板的水收集系數、以及三維旋轉航空發動機進口整流罩帽的水收集系數，并將數值計算結果與改進后的冰風洞運行相似準則結果進行對比，結果表明翼型的水收集系數和縮放后的翼型水收集系數吻合較好；張麗芬等[10]利用歐拉-拉格朗日法計算空氣-液滴兩相流，結合溢流水的溢流方向以及空氣速度矢量，整合出一套計算三維翼型表面積冰的方法，利用該方法對NACA0012平直翼和截面為GLC-305的后掠翼翼型進行非穩態、非結構網格的數值模擬，其結果與試驗值吻合較好。

目前對于結冰方面的研究主要集中在冰形預測上，對于結冰前后葉型氣動特性的對比分析相對較少。雖然有少數研究人員分析了結冰前后翼型、葉型周圍流場的變化，但缺少相關氣動參數的分析[11-12]。

本文通過商業軟件FENSAP-ICE模擬NACA0012翼型的結冰過程，對比NASA的結冰試驗數據，驗證結冰數值計算結果的正確性；對某型壓氣機進口導葉進行結冰計算，定量分析結冰前后葉片氣動參數的變化情況，總結結冰對葉型氣動特性的影響情況，以期為后續相關工作奠定基礎。

1 結冰數值模擬方法

本文采用定常描述計算結冰，假定空氣工質為理想氣體，忽略質量力。利用微小顆粒偏微分方程計算液滴速度，采用歐拉法計算液滴撞擊特性；利用Shallow-Water結冰模型，根據守恒定理對機翼表面的結冰量和結冰冰形進行求解。利用模塊化思想，使用時間多步長法，將總的結冰過程分成N步進行計算，在每個步長內分別進行繞流流場、液滴撞擊特性、熱力學結冰的計算，單步完成后重新劃分網格，繼續計算直至總的結冰過程結束。

1.1 微小顆粒偏微分方程

FENSAP-ICE解決了微小顆粒偏微分方程對于顆粒速度和水濃度的問題，因此可以計算獲得單次液滴噴射中的水濃度、液滴速度、液態水捕捉效率分布、撞擊模式和整個區域的撞擊極限，而不需要在噴射點處對液滴進行繁瑣的迭代，從而加快了液滴撞擊特性的計算效率。

1.2 Shallow-Water結冰模型

Shallow-Water結冰模型是一種基于表面水膜運動而建立的固體表面結冰模型，當液態水撞擊到固體表面時形成一層薄膜，在氣流的作用下，水膜會發生回流，加之熱力學的作用，水膜會發生凍結、蒸發和升華。模型通過確定固體表面每個節點上的水膜厚度，對結冰時和液態水回流時的傳熱和傳質問題進行計算。

Shallow-Water結冰模型假設水膜以液態水或積冰的形式附著在固體表面，在固體表面和法向上建立三維坐標系，并在該坐標系內構建水膜速度函數。同時，對水膜速度函數問題進行簡化，假設水膜速度沿固體表面法線方向呈線性變化，在粘性壁面處速度為0。對水膜上的速度進行平均，得到水膜厚度與水膜平均速度間的函數關系。最終通過水膜運動偏微分方程，即水膜的質量守恒方程和能量守恒方程計算固體表面的結冰量和結冰冰形。

2 NACA0012翼型結冰數值模擬

2.1 模型及邊界條件介紹

以結冰計算中常見的NACA0012翼型為計算模型，其翼型弦長0.533 4 m，在氣流攻角4°，空氣總壓101 325 Pa，空氣流速67 m/s，液態水含量1.0 g/m3，液滴顆粒直徑20 μm，結冰時間6 min的情況下，給定翼型壁面為固體無滑移壁面，分別對翼型在環境溫度為-2.22 ℃、-26.11 ℃時的結冰情況進行計算。其中，美國NASA的研究人員經過長期試驗和計算研究，對比冰形發現以1 min為間隔進行時間多步長法數值計算最為合理[13]，因此使用時間多步長法計算的時間間隔也選取1 min。

2.2 網格劃分

采用ICEM軟件對模型進行網格劃分，在翼型周圍采用O型結構化網格，計算模型翼型前緣網格分布情況如圖1所示，其中第一層網格高度0.003 mm，膨脹比1.1，保證網格y+值小于1。

圖1 NACA0012翼型計算網格

在空氣流速67 m/s，液態水含量1 g/m3，溫度-2.22 ℃，結冰時間6 min，攻角4°，液滴顆粒直徑20 μm的條件下，對計算模型進行網格無關性驗證，網格數量對翼型結冰量的影響規律分布如圖2所示，最終選擇20萬網格量進行計算。

圖2 網格數量與翼型結冰量關系

2.3 數值計算結果分析

2.3.1 冰形計算結果

環境溫度為-2.22 ℃、-26.11 ℃時翼型的結冰冰形計算結果如圖3所示，其中c為弦長。

(a) -2.22 ℃

從圖3可以看出：在其他條件一致的情況下，當環境溫度為較高的-2.22 ℃時，由于撞擊到翼型表面的液滴僅有一部分凍結為冰，另一部分受來流影響繼續向翼型的上方或后方移動最終凍結，這一部分溢流水在翼型前緣形成了向上翹起的角冰，因此翼型在環境溫度為-2.22 ℃時凝結成了典型的明冰；當環境溫度為較低的-26.11 ℃時，由于來流液滴撞擊到翼型表面后立即全部凍結，并無溢流水的形成，因此翼型在環境溫度為-26.11 ℃時，凍結成為外形與翼型型線較為一致，形狀較規則的毛冰。

2.3.2 計算結果誤差分析

為考察數值模擬方法的準確性，本文利用商業軟件FENSAP-ICE計算得到的翼型結冰冰形與現有的試驗冰形、結冰代碼預測冰形進行對比。試驗冰形來源于NASA Lewis研究中心的結冰風洞試驗臺[14-15]。Lewis研究中心利用自行研發的2D LEWICE/IBL代碼對NACA0012翼型結冰情況進行預測，計算數據豐富。為了精準地對翼型結冰計算冰形與試驗冰形進行定量分析，獲得結冰數值模擬方法的計算精度，引入三個無因次結冰特性量(無單位)對翼型結冰數值計算結果進行誤差分析，如表1～表2所示。對于NACA0012這類對稱翼型特征量[16]有三個。

表1 -2.22 ℃時翼型結冰特征量計算結果與誤差分析

表2 -26.11 ℃時翼型結冰特征量計算結果與誤差分析

無因次結冰面積ηs：

(1)

無因次結冰上極限ηLu：

(2)

無因次結冰下極限ηLd：

(3)

式中：A為翼型的面積，m2；c為翼型的弦長，m；S為翼型結冰面積，m2；Lu為結冰上極限，m；Ld為結冰下極限，m。

翼型結冰特征量示意圖如圖4所示。

圖4 翼型結冰特征量

從表1～表2可以看出：在環境溫度為-2.22 ℃和-26.11 ℃時，利用時間多步長法計算獲得的冰形在無因次結冰面積、無因次結冰上極限、無因次結冰下極限這三個結冰特征量的誤差整體小于由2D LEWICE/IBL代碼獲得的翼型結冰冰形誤差，主要是由于FENSAP-ICE將氣流粘性考慮在流場計算中，并利用Shallow-Water結冰模型進行冰形計算；結合翼型結冰冰形和結冰特征量誤差分析來看，利用FENSAP-ICE時間多步長數值模擬方法可以獲得與試驗冰形吻合情況較好的明冰、毛冰冰形，能夠清晰的反映出冰形的增長趨勢和變化情況，計算結果處于研究所允許的誤差范圍內。因此，本文采用的數值模擬方法獲得的明冰和毛冰冰形可用于后續研究。

3 壓氣機進口導葉結冰數值模擬

3.1 模型及邊界條件

以某型號壓氣機進口導葉為模型，不考慮葉片根部、冠部區域以及相鄰葉片間的影響，計算壓氣機進口導葉50%葉高處的結冰情況。其邊界條件為：氣流攻角0°，空氣總壓101 325 Pa，空氣流速100 m/s，液態水含量0.25 g/m3，液滴顆粒直徑10 μm，結冰時間60 s，多步長計算時間間隔20 s，葉型壁面為固體無滑移壁面；環境總溫分別為265 K和276 K。

3.2 網格劃分

使用ICEM軟件對模型進行網格劃分，在葉型周圍采用O型結構化網格，進口導葉葉型前緣、尾緣網格分布情況如圖5所示。

(a) 葉型前緣網格分布

在葉型計算網格中，第一層網格高度0.001 mm，膨脹比1.1。保證網格y+值小于1，在空氣流速100 m/s，液態水含量0.25 g/m3，溫度276 K，結冰時間60 s，攻角0°，液滴顆粒直徑10 μm的條件下，對計算模型進行網格無關性驗證，網格數量對葉型結冰量的影響規律分布如圖6所示，本文選擇20萬網格量進行計算。

圖6 網格數量與葉型結冰量關系

3.3 數值計算結果分析

3.3.1 冰形計算結果

環境溫度為265、276 K時葉型結冰冰形的計算結果如圖7所示。

(a) 265 K葉型結冰冰形

從圖7可以看出：當環境溫度為較低的265 K時，葉片前緣形成了與葉片形線較為一致的毛冰；當環境溫度為較高的276 K時，葉片前緣形成了向上翹起的明冰。這與NACA0012對稱翼型的結冰情況相同，對于具有非對稱結構的進口導葉來說，除了葉片前緣產生積冰外，在葉片壓力面上也發生了結冰現象。但壓力面上的積冰厚度不大，且冰形與葉片壓力面型線基本保持一致。

3.3.2 結冰前后流場分析

采用FENSAP-ICE計算葉型結冰，該軟件會自動生成由于結冰而發生位移后的葉型計算網格，如圖8所示。

(a) 265 K葉型前緣 (b) 265 K葉型尾緣

利用商業軟件Fluent，設置與初始流場計算相同的邊界條件，對結冰后葉型進行流場計算。結冰前后葉型的流場情況如圖9所示，可以看出：即使氣流以0°攻角流向葉型，但葉型的滯止區域仍發生在葉型的前緣偏下部，主要是由于氣流流過葉片吸力面的速度要大于氣流流過葉片壓力面的速度，因此在葉型前緣處生成了一個下大上小的壓力梯度，導致氣流以一定角度流向葉型，在結冰前后葉型的尾部上側都存在一個低速區，這個低速區內存在一個順時針的大面積分離渦和一個逆時針的小面積分離渦，其中順時針的分離渦主要是由于吸力面上發生的氣流分離造成的，而較小面積的逆時針分離渦是由于順時針分離渦與主流發生摻混而形成的；當環境溫度為265 K時，結冰后葉型僅在前緣上部形成了一片較小的氣流分離區，但并沒有影響葉片周圍流場的分布情況，結冰前后葉型的氣流分離都發生在葉片吸力面側，起始于弦長的50%處，毛冰的形成并沒有對葉片周圍流場產生較大影響；當環境溫度在276 K時，由于向上翹起的角冰，導致結冰后葉型的氣流分離加劇，分離點提前至葉片弦長的1%處，產生了一個面積約為原來2倍的低速分離區，低速區內的兩個分離渦面積也隨之增大，導致葉型的氣動性能急劇惡化。

(a) 265 K結冰前

結冰前后葉型出口處，即x=0.06 mm處氣流速度大小分布如圖10所示。當環境溫度為265 K時，結冰前后速度值變化不大。此時，低速區主要集中在-0.01～0.0 025 mm，這與圖9(a)，(b)所展現的馬赫云圖分布一致。但當環境溫度為276 K時，葉型出口處的最低速度值下降了約80%，低速區范圍擴大了40%。同時，當環境溫度為276 K時，結冰后葉型尾部低速區兩側的氣流速度明顯高于結冰前，可能是由于結冰導致的氣流分離加劇，渦流流速增加，摻混主流后造成的。此外，結冰后尾流低速區中部存在一個速度值先增大后減小的過程，而結冰前卻沒有這個過程，這主要是葉型前緣結冰，使得尾部分離區沿弦向擴大，分離區中兩個方向不同的分離渦在出口處共同作用而導致的氣流速度值增大。綜上，結冰會明顯增大葉型出口速度不均勻性分析。

(a) 265 K結冰前后葉型出口處速度分布

3.3.3 結冰前后總壓損失系數分析

葉型總壓損失系數ω的計算公式為

(4)

式中:PT1為葉型進口，即x=-0.02 mm處氣流總壓；PT2為葉型進口氣流動壓；PD1為葉型出口，即x=0.06 mm處氣流總壓。

不同環境溫度下葉型總壓損失系數分布如圖11所示。

(a) 265 K結冰前后葉型總壓損失系數分布

(b) 276 K結冰前后葉型總壓損失系數分布

從圖11(a)可以看出：環境溫度為265 K時，在y=-0.00 6 mm處總壓損失系數最大，整體來看，結冰后葉型的總壓損失系數有所增大，但增量極小。因此，葉型在265 K凍結形成的毛冰并不會額外產生較大的總壓損失。從圖11(b)可以看出：當環境溫度為276 K時，結冰后總壓損失系數明顯升高，在y=-0.006 mm處總壓損失系數達到最大值，較結冰前升高了約22%，隨后總壓損失系數呈現先下降，后上升再下降的趨勢，主要是由于出口處速度先上升后下降再上升而導致的。還可以看出：無論環境溫度高低，結冰后葉型壓力損失系數的增大主要集中在尾流分離區的上部，這表明，結冰后葉型的壓力損失大小主要由分離區上部的順時針分離渦影響。葉型前緣結冰上翹角度越大，葉型尾部分離區域面積越大，分離渦越強烈，導致結冰后葉型壓力損失越大。

4 結 論

(1) 明冰對于翼型和葉型的氣動影響均大于毛冰。

(2) 明冰會導致葉型吸力面上分離點提前，尾部分離區面積明顯增大，進而惡化葉型出口處速度分布，導致葉型總壓損失升高。而毛冰對葉型氣動特性影響很小。

(3) 葉型氣動特性的衰退主要受分離區中上分離渦影響，下分離渦則是由上分離渦和主流摻混而形成的，對葉型氣動特性影響不大。

(4) 采用商業軟件FENSAP-ICE可以高效、準確地對明冰、毛冰這兩種典型的翼型、葉型結冰情況進行預測。