姜 寒，朱春玲

(南京航空航天大學(xué)航空學(xué)院，南京 210016)

飛機(jī)機(jī)翼作為提供升力的重要部件，在飛行過程中起著至關(guān)重要的作用[1]。當(dāng)飛機(jī)在空中遇到結(jié)冰氣象條件時，機(jī)翼表面會發(fā)生結(jié)冰現(xiàn)象，從而改變機(jī)翼的氣動外形，導(dǎo)致升力下降阻力升高，影響飛行安全。因此對機(jī)翼結(jié)冰的研究就尤為重要，目前對機(jī)翼結(jié)冰研究的方法主要有結(jié)冰風(fēng)洞實驗和數(shù)值模擬兩種。

結(jié)冰風(fēng)洞實驗主要通過將模型固定在結(jié)冰風(fēng)洞內(nèi)，通過風(fēng)機(jī)和噴霧在制冷條件下模擬結(jié)冰氣象條件下的飛行過程，觀察結(jié)冰情況。

數(shù)值模擬利用計算機(jī)實現(xiàn)不同條件下的流場水滴場以及結(jié)冰的計算，因其應(yīng)用的自由度和靈活性成為結(jié)冰研究中的主要方式。國外針對結(jié)冰計算開發(fā)了多款軟件，例如 FENSAP-ICE[2]、LEWICE[3]和 ONERA[4]等。其中一些軟件對結(jié)冰的模擬主要基于 Messinger 模型[5]。該模型通過求解表面積冰控制體的質(zhì)量守恒及能量守恒方程，得到控制體的平衡溫度，進(jìn)而求得結(jié)冰量。Myers等[6]提出了能在結(jié)冰增長過程實現(xiàn)水膜、冰層之間傳熱傳質(zhì)的計算模型。中國對于結(jié)冰計算的研究工作也在發(fā)展中，張強(qiáng)等[7]基于歐拉兩相流理論進(jìn)行了三維機(jī)翼表面的霜冰模擬；曹廣州等[8]根據(jù)冰層表面水膜流動特點(diǎn)建立了能模擬水膜流動的三維積冰數(shù)學(xué)模型；常士楠等[9]利用歐拉法求解空氣水滴兩相流流場并改進(jìn)了 Messinger 模型；雷夢龍等[10]通過修改結(jié)冰類型判斷方法，改進(jìn)了 Myers 模型；王翊等[11]通過數(shù)值模擬的方法研究了二維翼型冰型生長過程；王強(qiáng)等[12]對二維風(fēng)力機(jī)翼型進(jìn)行了樹脂及實驗研究?？梢钥吹剑陙碇袊鴮W(xué)者對于機(jī)翼結(jié)冰的數(shù)值及實驗研究越來越多，但大多數(shù)仍集中于二維翼型，對后掠翼結(jié)冰研究較少。

在建立數(shù)值計算模型過程中，現(xiàn)考慮水滴撞擊在機(jī)翼表面先形成一層水膜，基于這層水膜進(jìn)行結(jié)冰增長的計算，實現(xiàn)三維后掠翼的結(jié)冰模擬。設(shè)計后掠翼三維模型，在結(jié)冰風(fēng)洞進(jìn)行實驗，與計算結(jié)果進(jìn)行對比，驗證計算模型的可行性，觀察后掠翼積冰外形，得到后掠翼沿展向結(jié)冰的規(guī)律。

1 后掠翼結(jié)冰數(shù)值模擬

1.1 三維結(jié)冰模型的建立

首先，對模型進(jìn)行如下假設(shè)。

(1)撞擊在機(jī)翼表面的水滴會先在重力、剪切力及壓力梯度的作用下向下游流動形成一層穩(wěn)態(tài)的初始水膜，后結(jié)冰。

(2)結(jié)冰時忽略物面與氣流間的輻射換熱。

(3)未凍結(jié)的水在向下游流動時會形成連續(xù)的水膜且不發(fā)生破裂。

基于以上假設(shè)，建立初始水膜流動模型。建立貼體坐標(biāo)系(ξ，η，ζ)，其中ξ、η為曲面上的參考坐標(biāo)，ζ垂直于ξ、η。

機(jī)翼表面初始水膜流動過程中，通過控制體的質(zhì)量包括：水滴撞擊項mimp，水蒸發(fā)項mevp，通過控制體前后左右表面的流動項。對控制體應(yīng)用質(zhì)量守恒定理得到水膜連續(xù)性方程為

(1)

式(1)中：u、v分別為ξ、η方向的速度。

mimp=βV∞CLW

(2)

式(2)中：β為水收集系數(shù)；V∞為來流速度，m/s；CLW為液態(tài)水含量，kg/m3。

(3)

式(3)中：hc為對流換熱系數(shù)，W/(m2·K)；ρa(bǔ)、cp、Rv分別為水蒸氣密度、定壓比熱容及氣體常數(shù)；esur、e∞分別為表面和大氣飽和水蒸氣分壓力，Pa；Tsur、T∞分別為表面和大氣溫度，K；φ為大氣相對濕度。

通過量綱分析法簡化水膜動量方程為

(4)

式(4)中：P為空氣壓力；μw為水動力黏性系數(shù)；ρw為水密度；gξ、gη分別為重力在ξ、η方向的分量。

邊界采用無滑移邊界條件，即

(5)

式(5)中：τa為水膜表面剪切應(yīng)力；Hw為水膜厚度。

基于上述初始水膜模型，建立結(jié)冰增長模型。對每一個控制體建立水的質(zhì)量守恒方程為

mini+mimp+min-mout-mevp-mice=0

(6)

式(6)中：mini為初始水膜項，mini=ρwHw；min、mout分別為流入、流出項；mice為結(jié)冰項。

建立水的能量守恒方程為

qimp+qin+qcnd-qevp-qice-qcnv-qout=0

(7)

式(7)中：qimp為撞擊項；qin、qout分別為流入、流出項；qcnd為導(dǎo)熱項；qice為結(jié)冰項；qcnv為對流換熱項；qevp為蒸發(fā)項。

各項熱流計算式為

(8)

qin=min[ciTm+LF+cw(Tin-Tm)]

(9)

qout=mout[ciTm+LF+cw(Tsur-Tm)]

(10)

qice=miceciTsur

(11)

qcnv=hc(Tsur-T∞)

(12)

qevp=mevp{ciTm+LF+[frzci+(1-frz)cw]×

(Tsur-Tm)+LE}

(13)

式中：Tsur、T∞分別為表面溫度和來流溫度；Tw、Tin、Tm分別為水膜溫度、入流水溫及相變平衡溫度，K；LF、LE、LS分別為冰的融解潛熱、水的蒸發(fā)潛熱、冰升華潛熱，J/kg；cw、ci分別為水和冰的比熱；frz為引入的凍結(jié)系數(shù)，表征凍結(jié)水的質(zhì)量與進(jìn)入控制體中水質(zhì)量的比值，計算公式為

frz=mice/(mini+mimp+min)

(14)

1.2 三維結(jié)冰模型的求解

求解初始水膜模型時，先利用邊界條件式(5)求解水膜動量方程(4)，可得到水膜沿法線方向的速度分布為

(15)

(16)

分別對式(15)、式(16)在η方向積分然后取平均，即可得到水膜沿翼面流動時分別沿兩個方向ξ、η的平均速度為

(17)

(18)

將式(17)、式(18)代入式(1)，可得到關(guān)于水膜厚度Hw的偏微分方程為

(19)

(1)將所有控制體積中水的流入項賦初值為0。

(2)聯(lián)立水的質(zhì)量和能量守恒方程，代入式(14)，得到凍結(jié)系數(shù)的表達(dá)式，從而計算凍結(jié)系數(shù)frz，可得

frz={qclt+qin+qcnd-qevp-qcnv-(mini+mclt+min-mevp)[ciTm+LF+cw(Tsur-Tm)]}/{(mini+mclt+min)[(ci-cw)(Tsur-Tm)-LF]}

(20)

(3)用式(21)計算流出項mout，并根據(jù)水膜模型中計算出水流動速度分量按比例分配各方向的流出分量。

mout=(1-frz)(mini+mclt+min)-mevp

(21)

(4)根據(jù)流出分量重新計算各控制體的流入項并計算流出水的殘差值e，公式為

e=max[mout(i)-mout(i-1)]

(22)

式(22)中：i為當(dāng)前時間步；i-1為上一時間步。

(5)判斷流出水項殘差e是否滿足收斂條件，若不收斂，重復(fù)上述步驟至收斂。

2 后掠翼結(jié)冰試驗

2.1 實驗?zāi)Ｐ团c測量方法

設(shè)計后掠角為 30°的 NACA0012 機(jī)翼模型，使其能與結(jié)冰風(fēng)洞支架實現(xiàn)連接，由于試驗段寬300 mm，將模型展長定為 340 mm，使其在安裝進(jìn)試驗段后與兩邊玻璃有一點(diǎn)空隙，NACA0012截面弦長定為 300 mm。利用 3D 打印技術(shù)加工該模型。進(jìn)行實驗時，將模型通過支架固定在結(jié)冰風(fēng)洞試驗段中，結(jié)冰試驗結(jié)束后，將帶冰的機(jī)翼模型連同支架一同從試驗段中取出放在低溫環(huán)境中以保證冰型不被破壞。結(jié)冰試驗臺如圖1所示。

圖1 結(jié)冰試驗臺

為準(zhǔn)確測得實驗所得結(jié)冰冰型，需要利用三維掃描儀先對未結(jié)冰的干凈帶支架翼型進(jìn)行掃描，實驗后在低溫環(huán)境下對結(jié)冰帶支架翼型進(jìn)行掃描，將二者根據(jù)結(jié)冰實驗前后完全沒發(fā)生變化的支架上3個不同方向特征面進(jìn)行重合處理，以得到二者相差部分即為機(jī)翼表面結(jié)冰的外形。

2.2 后掠翼結(jié)冰試驗與結(jié)果

為研究機(jī)翼表面存在溢流水時的結(jié)冰情況，將來流溫度控制在-8 ℃，工況如表1所示。

表1 結(jié)冰工況

結(jié)冰結(jié)果如圖2所示，4個結(jié)冰截面厚度的掃描儀還原模型如圖3所示。

圖2 結(jié)冰試驗結(jié)果

圖3 掃描儀還原冰型

在還原模型上截取4個截面如圖4所示，以便后期與結(jié)冰計算結(jié)果進(jìn)行對比。圖4(a)為截面z=60 mm處的結(jié)冰冰型，圖4(b)為截面z=110 mm 處的結(jié)冰冰型，圖4(c)為截面z=190 mm 處的結(jié)冰冰型，圖4(d)為截面z=240 mm 處的結(jié)冰冰型。由圖4可知：在攻角為 4°時，下翼面積冰較多且積冰范圍較大，冰型輪廓較為光滑，無明顯冰角凸起。上翼面積冰較少且積冰范圍較小，有冰角，符合明冰結(jié)冰特征。由于結(jié)冰溫度選取了溫度較低的結(jié)明冰工況，因此冰角并沒有典型明冰工況實驗結(jié)果的冰角突出。圖5為由翼根至翼梢各截面的結(jié)冰冰型對比?？梢钥闯?，沿展向方向由翼根至翼梢結(jié)冰冰型輪廓相似但結(jié)冰量逐漸增加，這是由于后掠角帶來的氣流展向流動，使水滴在撞擊在機(jī)翼上形成的水膜除了弦向的流動外，還會沿展向流動，導(dǎo)致翼梢積冰量大于翼根。

圖4 二維冰型截面圖

圖5 試驗結(jié)果截面對比圖

3 算例驗證結(jié)果

為驗證本文三維結(jié)冰模型的正確性，選擇實驗采用的后掠角為30°的NACA0012機(jī)翼作為計算模型，為與實驗結(jié)果進(jìn)行對比，對弦長為1 m的標(biāo)準(zhǔn)計算模型進(jìn)行縮比至實驗?zāi)Ｐ拖嗤叽?，弦長為0.3 m。其幾何外形如圖6所示。

圖6 NACA0012 30°后掠機(jī)翼幾何外形

結(jié)冰計算狀態(tài)采用與實驗相同的工況，如表1所示。

圖7為結(jié)冰計算結(jié)果的三維視圖，圖中紅色網(wǎng)格為干凈翼面。沿翼型展向取4個截面，與FENSAP-ICE軟件數(shù)值模擬結(jié)果與實驗結(jié)果進(jìn)行對比，選取的4個截面如圖7所示。

圖7 三維冰型圖

圖8為截面z分別為60、110、190、240 mm處的結(jié)冰冰型對比。從圖8可以看出：本文模型結(jié)果與FENSAP-ICE結(jié)冰輪廓基本吻合。與實驗結(jié)果對比可以發(fā)現(xiàn)，在靠近翼根的截面z=60 mm 處，本文模型求得的結(jié)冰量略小于實驗結(jié)果，冰型輪廓基本符合，冰角形狀不明顯。在翼段中部的截面z=110 mm 和截面z=190 mm 處，本文模型計算結(jié)果與實驗結(jié)果較為接近，對于冰角形狀的計算略有不足。而在靠近翼梢的截面z=240 mm 處，本文模型求解的結(jié)冰量大于實驗結(jié)果，冰角位置較為相似。

圖8 截面冰型對比圖

圖9為數(shù)值計算結(jié)果各截面的冰型對比，從圖9中可以看出，同實驗結(jié)果相同，由翼根至翼梢，結(jié)冰量逐漸增長，但模型對冰角形狀的計算不是很準(zhǔn)確，原因可能是在計算結(jié)冰增長過程中，采用準(zhǔn)穩(wěn)態(tài)方法計算，雖考慮結(jié)冰初始水膜流動，但對其在結(jié)冰增長過程中由于表面形狀變化導(dǎo)致的流動沒有考慮，導(dǎo)致了結(jié)冰特征外形如冰角形狀的不準(zhǔn)確。

圖9 計算結(jié)果截面對比圖

4 結(jié)論

考慮水滴撞擊在機(jī)翼表面形成的初始水膜，求解結(jié)冰增長過程的質(zhì)量與能量守恒方程得到能夠計算后掠翼結(jié)冰的數(shù)值方法。在結(jié)冰風(fēng)洞中進(jìn)行 30°后掠角的 NACA0012 三維模型結(jié)冰實驗，利用三維掃描儀還原冰型。將計算結(jié)果與FENSAP-ICE 結(jié)果以及實驗結(jié)果進(jìn)行對比，得出以下結(jié)論。

(1)本文模型考慮了明冰工況下水滴撞擊翼面形成的初始水膜流動，并將算得參數(shù)應(yīng)用在結(jié)冰增長過程中，與考慮了結(jié)冰過程中水膜流動的FENSAP-ICE 軟件計算結(jié)果符合良好。

(2)結(jié)冰風(fēng)洞實驗結(jié)果符合明冰工況結(jié)冰冰型，說明實驗結(jié)果可用來驗證計算。將計算結(jié)果與之對比可以發(fā)現(xiàn)基本吻合，在冰角處形狀略有差異，證明本文計算方法可行。

(3)后掠翼模型在結(jié)明冰時可以看到由翼根至翼梢結(jié)冰厚度逐漸變大，說明后掠翼結(jié)冰過程中水沿展向的流動會對結(jié)冰外形產(chǎn)生明顯的影響。