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

(1. 哈爾濱工程大學(xué) 動(dòng)力與能源工程學(xué)院， 哈爾濱 150001)

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

0 引 言

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

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

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

目前對(duì)于結(jié)冰方面的研究主要集中在冰形預(yù)測(cè)上，對(duì)于結(jié)冰前后葉型氣動(dòng)特性的對(duì)比分析相對(duì)較少。雖然有少數(shù)研究人員分析了結(jié)冰前后翼型、葉型周圍流場(chǎng)的變化，但缺少相關(guān)氣動(dòng)參數(shù)的分析[11-12]。

本文通過商業(yè)軟件FENSAP-ICE模擬NACA0012翼型的結(jié)冰過程，對(duì)比NASA的結(jié)冰試驗(yàn)數(shù)據(jù)，驗(yàn)證結(jié)冰數(shù)值計(jì)算結(jié)果的正確性；對(duì)某型壓氣機(jī)進(jìn)口導(dǎo)葉進(jìn)行結(jié)冰計(jì)算，定量分析結(jié)冰前后葉片氣動(dòng)參數(shù)的變化情況，總結(jié)結(jié)冰對(duì)葉型氣動(dòng)特性的影響情況，以期為后續(xù)相關(guān)工作奠定基礎(chǔ)。

1 結(jié)冰數(shù)值模擬方法

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

1.1 微小顆粒偏微分方程

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

1.2 Shallow-Water結(jié)冰模型

Shallow-Water結(jié)冰模型是一種基于表面水膜運(yùn)動(dòng)而建立的固體表面結(jié)冰模型，當(dāng)液態(tài)水撞擊到固體表面時(shí)形成一層薄膜，在氣流的作用下，水膜會(huì)發(fā)生回流，加之熱力學(xué)的作用，水膜會(huì)發(fā)生凍結(jié)、蒸發(fā)和升華。模型通過確定固體表面每個(gè)節(jié)點(diǎn)上的水膜厚度，對(duì)結(jié)冰時(shí)和液態(tài)水回流時(shí)的傳熱和傳質(zhì)問題進(jìn)行計(jì)算。

Shallow-Water結(jié)冰模型假設(shè)水膜以液態(tài)水或積冰的形式附著在固體表面，在固體表面和法向上建立三維坐標(biāo)系，并在該坐標(biāo)系內(nèi)構(gòu)建水膜速度函數(shù)。同時(shí)，對(duì)水膜速度函數(shù)問題進(jìn)行簡(jiǎn)化，假設(shè)水膜速度沿固體表面法線方向呈線性變化，在粘性壁面處速度為0。對(duì)水膜上的速度進(jìn)行平均，得到水膜厚度與水膜平均速度間的函數(shù)關(guān)系。最終通過水膜運(yùn)動(dòng)偏微分方程，即水膜的質(zhì)量守恒方程和能量守恒方程計(jì)算固體表面的結(jié)冰量和結(jié)冰冰形。

2 NACA0012翼型結(jié)冰數(shù)值模擬

2.1 模型及邊界條件介紹

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

2.2 網(wǎng)格劃分

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

圖1 NACA0012翼型計(jì)算網(wǎng)格

在空氣流速67 m/s，液態(tài)水含量1 g/m3，溫度-2.22 ℃，結(jié)冰時(shí)間6 min，攻角4°，液滴顆粒直徑20 μm的條件下，對(duì)計(jì)算模型進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證，網(wǎng)格數(shù)量對(duì)翼型結(jié)冰量的影響規(guī)律分布如圖2所示，最終選擇20萬網(wǎng)格量進(jìn)行計(jì)算。

圖2 網(wǎng)格數(shù)量與翼型結(jié)冰量關(guān)系

2.3 數(shù)值計(jì)算結(jié)果分析

2.3.1 冰形計(jì)算結(jié)果

環(huán)境溫度為-2.22 ℃、-26.11 ℃時(shí)翼型的結(jié)冰冰形計(jì)算結(jié)果如圖3所示，其中c為弦長(zhǎng)。

(a) -2.22 ℃

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

2.3.2 計(jì)算結(jié)果誤差分析

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

表1 -2.22 ℃時(shí)翼型結(jié)冰特征量計(jì)算結(jié)果與誤差分析

表2 -26.11 ℃時(shí)翼型結(jié)冰特征量計(jì)算結(jié)果與誤差分析

無因次結(jié)冰面積ηs：

(1)

無因次結(jié)冰上極限ηLu：

(2)

無因次結(jié)冰下極限ηLd：

(3)

式中：A為翼型的面積，m2；c為翼型的弦長(zhǎng)，m；S為翼型結(jié)冰面積，m2；Lu為結(jié)冰上極限，m；Ld為結(jié)冰下極限，m。

翼型結(jié)冰特征量示意圖如圖4所示。

圖4 翼型結(jié)冰特征量

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

3 壓氣機(jī)進(jìn)口導(dǎo)葉結(jié)冰數(shù)值模擬

3.1 模型及邊界條件

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

3.2 網(wǎng)格劃分

使用ICEM軟件對(duì)模型進(jìn)行網(wǎng)格劃分，在葉型周圍采用O型結(jié)構(gòu)化網(wǎng)格，進(jìn)口導(dǎo)葉葉型前緣、尾緣網(wǎng)格分布情況如圖5所示。

(a) 葉型前緣網(wǎng)格分布

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

圖6 網(wǎng)格數(shù)量與葉型結(jié)冰量關(guān)系

3.3 數(shù)值計(jì)算結(jié)果分析

3.3.1 冰形計(jì)算結(jié)果

環(huán)境溫度為265、276 K時(shí)葉型結(jié)冰冰形的計(jì)算結(jié)果如圖7所示。

(a) 265 K葉型結(jié)冰冰形

從圖7可以看出：當(dāng)環(huán)境溫度為較低的265 K時(shí)，葉片前緣形成了與葉片形線較為一致的毛冰；當(dāng)環(huán)境溫度為較高的276 K時(shí)，葉片前緣形成了向上翹起的明冰。這與NACA0012對(duì)稱翼型的結(jié)冰情況相同，對(duì)于具有非對(duì)稱結(jié)構(gòu)的進(jìn)口導(dǎo)葉來說，除了葉片前緣產(chǎn)生積冰外，在葉片壓力面上也發(fā)生了結(jié)冰現(xiàn)象。但壓力面上的積冰厚度不大，且冰形與葉片壓力面型線基本保持一致。

3.3.2 結(jié)冰前后流場(chǎng)分析

采用FENSAP-ICE計(jì)算葉型結(jié)冰，該軟件會(huì)自動(dòng)生成由于結(jié)冰而發(fā)生位移后的葉型計(jì)算網(wǎng)格，如圖8所示。

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

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

(a) 265 K結(jié)冰前

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

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

3.3.3 結(jié)冰前后總壓損失系數(shù)分析

葉型總壓損失系數(shù)ω的計(jì)算公式為

(4)

式中:PT1為葉型進(jìn)口，即x=-0.02 mm處氣流總壓；PT2為葉型進(jìn)口氣流動(dòng)壓；PD1為葉型出口，即x=0.06 mm處氣流總壓。

不同環(huán)境溫度下葉型總壓損失系數(shù)分布如圖11所示。

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

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

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

4 結(jié) 論

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

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

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

(4) 采用商業(yè)軟件FENSAP-ICE可以高效、準(zhǔn)確地對(duì)明冰、毛冰這兩種典型的翼型、葉型結(jié)冰情況進(jìn)行預(yù)測(cè)。