郭琦,申曉斌,2,*,林貴平,張世娟

(1.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院,北京 100191;2.中國空氣動(dòng)力研究與發(fā)展中心 結(jié)冰與防除冰重點(diǎn)實(shí)驗(yàn)室,綿陽 621000; 3.武漢航空儀表有限責(zé)任公司,武漢 430074)

飛機(jī)在結(jié)冰氣象條件下飛行時(shí),其迎風(fēng)面上會(huì)積聚過冷水滴,發(fā)生冰的聚積,該現(xiàn)象為飛機(jī)結(jié)冰[1]。飛機(jī)積冰對(duì)飛機(jī)性能有重要影響,會(huì)造成升力和失速迎角減小,阻力增加[2],甚至導(dǎo)致嚴(yán)重的飛行事故[3]。隨著航空科學(xué)技術(shù)的快速發(fā)展,航空飛行器得到廣泛運(yùn)用,由于飛機(jī)結(jié)冰引發(fā)的飛行安全問題愈發(fā)凸顯。目前對(duì)飛機(jī)結(jié)冰的研究有結(jié)冰氣象條件的探測(cè)、過冷水滴的撞擊特性、結(jié)冰熱力學(xué)模型與冰形預(yù)測(cè)、結(jié)冰對(duì)飛行性能的影響等,本文針對(duì)飛機(jī)旋轉(zhuǎn)表面結(jié)冰進(jìn)行數(shù)值模擬研究。

目前主要通過實(shí)驗(yàn)研究與數(shù)值仿真分析,對(duì)飛機(jī)結(jié)冰展開研究,其中數(shù)值方法具有經(jīng)濟(jì)性好、執(zhí)行周期短等優(yōu)點(diǎn),并且由于計(jì)算流體力學(xué)技術(shù)和計(jì)算機(jī)技術(shù)的不斷發(fā)展與成熟,數(shù)值仿真在飛機(jī)結(jié)冰研究中的作用越發(fā)重要。在國內(nèi)外學(xué)者的共同努力下,飛機(jī)結(jié)冰數(shù)值模擬已取得了相當(dāng)大的進(jìn)展,結(jié)冰數(shù)值研究的內(nèi)容由霜冰預(yù)測(cè)發(fā)展到了明冰的模擬,由二維計(jì)算逐步向三維仿真推進(jìn),由固定部件表面向旋轉(zhuǎn)部件表面發(fā)展,三維旋轉(zhuǎn)表面復(fù)雜流動(dòng)結(jié)冰機(jī)理與模擬成為了目前研究的焦點(diǎn)。

傳統(tǒng)的結(jié)冰計(jì)算一般采用Messinger[4]結(jié)冰熱力學(xué)模型,該模型假設(shè)某個(gè)控制體中的未凝結(jié)液態(tài)水會(huì)全部流出到其他控制體中,沒有考慮到控制體自身存在水膜的情況。美國國家航空航天局(NASA)于20世紀(jì)80年代研發(fā)成功的LEWICE軟件[5]就是基于該結(jié)冰模型開發(fā)的。之后在Messinger模型的基礎(chǔ)上,發(fā)展出了Myers[6]模型和Shallow-Water[7]結(jié)冰熱力學(xué)模型,兩者都考慮了控制體表面水膜存在的情況,不同的是,Myers模型還考慮了冰層和水膜中的溫度梯度分布。加拿大的FENSAP-ICE軟件[7],是利用Shallow-Water模型對(duì)結(jié)冰表面的冰形進(jìn)行求解。基于這些結(jié)冰模型,國內(nèi)外學(xué)者對(duì)飛機(jī)結(jié)冰展開了一系列模擬仿真研究。易賢等[8]發(fā)展了一種采用迭代求解的三維結(jié)冰熱力學(xué)模型,該模型基于傳統(tǒng)的Messinger二維模型,引入了結(jié)冰表面溢流水的流動(dòng)。申曉斌等[9]基于Messinger結(jié)冰模型,建立了一種適用于三維結(jié)冰表面的液態(tài)水溢流模型,并完成了對(duì)三維發(fā)動(dòng)機(jī)結(jié)冰冰形的計(jì)算。曹廣州等[10]采用Myers結(jié)冰模型,對(duì)三維機(jī)翼模型的結(jié)冰增長和表面液態(tài)水流動(dòng)進(jìn)行了仿真求解,且能夠在計(jì)算中自動(dòng)辨別過冷水滴在結(jié)冰表面發(fā)生碰撞后的結(jié)冰冰形。Cao和Huang[11]在Myers模型的基礎(chǔ)上,建立了三維結(jié)冰模擬模型,計(jì)算了后掠機(jī)翼的冰形。雷夢(mèng)龍等[12]基于Myers結(jié)冰熱力學(xué)模型,對(duì)其結(jié)冰類型判別方法進(jìn)行改進(jìn),采用離散算法對(duì)表面液態(tài)水溢流和冰層生長過程進(jìn)行計(jì)算,模擬了機(jī)翼表面的結(jié)冰過程。

同時(shí),隨著飛機(jī)結(jié)冰模擬技術(shù)的發(fā)展,各國學(xué)者針對(duì)飛機(jī)旋轉(zhuǎn)部件結(jié)冰仿真計(jì)算也進(jìn)行了多項(xiàng)研究。Bain等[13]用LEWICE3D軟件進(jìn)行旋轉(zhuǎn)槳葉的結(jié)冰計(jì)算,將旋轉(zhuǎn)結(jié)冰過程分為2個(gè)步長進(jìn)行計(jì)算,沒有考慮旋轉(zhuǎn)作用對(duì)表面水膜的作用力及對(duì)流動(dòng)特性的影響。Reid等[14-15]基于Shallow-Water溢流模型,利用旋轉(zhuǎn)坐標(biāo)系的方法對(duì)旋轉(zhuǎn)槳葉的結(jié)冰進(jìn)行了模擬。Switchenko等[16]針對(duì)旋轉(zhuǎn)槳葉采用二維截面的方法簡化計(jì)算,通過二維多步法模擬特征截面處的結(jié)冰情況。Wang與Zhu[17]在旋轉(zhuǎn)坐標(biāo)系下計(jì)算了流場(chǎng)與水滴撞擊特性,考慮離心力對(duì)溢流水膜分配的影響,進(jìn)行了結(jié)冰分析。Dong等[18]采用歐拉法求解了發(fā)動(dòng)機(jī)處的水滴撞擊特性,根據(jù)結(jié)冰熱力學(xué)平衡模型得到了旋轉(zhuǎn)帽罩的結(jié)冰冰形。趙秋月等[19]基于拉格朗日法,對(duì)旋轉(zhuǎn)部件表面的水收集系數(shù)計(jì)算方法進(jìn)行研究,將其用于發(fā)動(dòng)機(jī)旋轉(zhuǎn)整流帽罩的水滴撞擊特性計(jì)算,發(fā)現(xiàn)整流罩表面的局部水收集系數(shù)會(huì)隨著轉(zhuǎn)速和錐體錐角的增大,沿表面外形線下降加快,撞擊區(qū)域減小。吳孟龍等[20]基于歐拉法對(duì)發(fā)動(dòng)機(jī)旋轉(zhuǎn)帽罩三維表面的水滴撞擊特性進(jìn)行模擬,研究帽罩表面的水滴撞擊特性受旋轉(zhuǎn)速度的影響,發(fā)現(xiàn)旋轉(zhuǎn)速度對(duì)帽罩表面的水滴撞擊特性影響較小。

現(xiàn)在使用的結(jié)冰熱力學(xué)模型幾乎都是在固定部件表面建立的,一般認(rèn)為空氣剪切力為水膜的驅(qū)動(dòng)力,用其來確定水膜的流動(dòng)特性,而對(duì)于旋轉(zhuǎn)部件表面的受力分析及旋轉(zhuǎn)作用對(duì)水膜運(yùn)動(dòng)的影響分析較少,很少考慮離心力和科里奧利力的方向與大小對(duì)結(jié)冰效果的作用,旋轉(zhuǎn)結(jié)冰過程中的運(yùn)動(dòng)力學(xué)與傳熱學(xué)機(jī)理研究尚不全面。本文采用Shallow-Water結(jié)冰模型,對(duì)飛機(jī)旋轉(zhuǎn)表面的水膜流動(dòng)和結(jié)冰冰形進(jìn)行模擬,并分析轉(zhuǎn)速、水滴直徑和液態(tài)水含量等因素對(duì)旋轉(zhuǎn)表面結(jié)冰的影響,為飛機(jī)旋轉(zhuǎn)部件的結(jié)冰研究提供參考。

1 數(shù)學(xué)模型

為了研究旋轉(zhuǎn)表面的水膜流動(dòng)和結(jié)冰冰形,基于Shallow-Water結(jié)冰熱力學(xué)模型[7],分析結(jié)冰表面水膜的流動(dòng)與傳熱過程,建立動(dòng)量守恒方程、質(zhì)量守恒方程和能量守恒方程,進(jìn)而計(jì)算旋轉(zhuǎn)坐標(biāo)系下的非穩(wěn)態(tài)結(jié)冰過程與冰形。

1.1 動(dòng)量守恒方程

旋轉(zhuǎn)表面的水膜主要受到空氣剪切力、壓力梯度、重力、旋轉(zhuǎn)導(dǎo)致的慣性力等因素的影響,選取表面微元控制體進(jìn)行受力分析,如圖1所示。用S表示曲面的表面,S′表示S在法向量n方向的投影,則S和S′之間的距離就是水膜的厚度h。

圖1 旋轉(zhuǎn)表面水膜受力分析示意圖Fig.1 Schematic diagram of force analysis of water film on rotating surface

由于水膜在結(jié)冰表面流動(dòng)速度較小,且不考慮水滴撞擊的影響,對(duì)水膜作常物性處理,可假設(shè)結(jié)冰表面的水膜流動(dòng)為不可壓層流流動(dòng)[10],基于不可壓流體的雷諾平均Navier-Stokes方程對(duì)水膜溢流過程進(jìn)行建模,在旋轉(zhuǎn)坐標(biāo)系下其動(dòng)量方程表述為

式中:ρw為水的密度;v為水膜相對(duì)速度;p為水膜內(nèi)部壓力;I為單位張量;μw為水的動(dòng)力黏度;y為表面法向坐標(biāo);g為重力加速度;ω為轉(zhuǎn)速;r為當(dāng)?shù)厥笍健Ｆ渲?方程右邊第2項(xiàng)為慣性力項(xiàng),包括了科里奧利慣性力與牽連慣性力。

考慮結(jié)冰表面的水膜厚度很小,其主要受空氣剪切力和旋轉(zhuǎn)導(dǎo)致的慣性力的作用,忽略非穩(wěn)態(tài)項(xiàng)、壓力梯度項(xiàng)和重力項(xiàng)的影響[21],則動(dòng)量方程可簡化為

水膜沿表面的流動(dòng)方向?yàn)閤=(x1,x2)。垂直方向水膜流動(dòng)的邊界分別是水膜-固體表面界面和空氣-水膜界面。

水膜-表面界面在旋轉(zhuǎn)參考系下的條件為

空氣-水膜界面處,由于表面水膜很薄,水膜的速度梯度由空氣剪切力占主導(dǎo)作用,邊界條件為

式中:τ為水膜所受的空氣剪切力。

根據(jù)動(dòng)量方程(2)和邊界條件(3)、(4),得到旋轉(zhuǎn)坐標(biāo)系下,水膜速度沿y方向的分布:

水膜和壁面交界處,速度為0;水膜頂端,速度達(dá)到最大值。在法向進(jìn)行積分平均化,得到水膜的相對(duì)速度為

可見,水膜速度是由慣性力加速度和剪切應(yīng)力共同決定的,而慣性力加速度與該速度本身有關(guān),存在耦合關(guān)系,在求解過程中需要通過耦合迭代求解的方式予以處理。

1.2 質(zhì)量守恒方程

圖2為所有流入與流出控制體的質(zhì)量,包括:控制體上表面水滴的撞擊和蒸發(fā)或升華,下表面冰的生成,同時(shí)還有水的流入流出,以及控制體自身水膜厚度的變化[10]。

根據(jù)圖2可得質(zhì)量守恒方程為

圖2 控制體質(zhì)量守恒示意圖Fig.2 Schematic diagram of mass conservation of control volume

式(7)可表示為

式中:v1和v2分別為水膜平均速度v在x1和x2方向分量的大小。

由于非穩(wěn)態(tài)計(jì)算中,是對(duì)每一個(gè)時(shí)間步長進(jìn)行分別計(jì)算,需要將結(jié)冰時(shí)間分成多個(gè)時(shí)間步長,每一個(gè)時(shí)間步長作為一個(gè)穩(wěn)態(tài)。因此,對(duì)質(zhì)量守恒方程進(jìn)行離散求解,可表示為

式中:hold為上一時(shí)間步長控制體內(nèi)的水膜厚度;A為控制體底面面積;Δt為時(shí)間步長;nn為流入流出面的單位外法向量;Δx1和Δx2分別為當(dāng)前控制體的底面邊長。

式(9)中,等號(hào)左邊第1項(xiàng)對(duì)應(yīng)的是控制體內(nèi)水膜質(zhì)量隨時(shí)間的變化,第2和第3項(xiàng)表示控制體內(nèi)液態(tài)水在x1和x2方向所有的流入流出。從而,可將每個(gè)時(shí)間步長結(jié)冰計(jì)算中的質(zhì)量守恒方程表示為

式中:m和m分別為上一時(shí)間步長控制體內(nèi)存在的水膜和冰層的質(zhì)量;mw和mice分別為當(dāng)前時(shí)間步長控制體內(nèi)存在的水膜和冰層的質(zhì)量;min和mout分別為當(dāng)前時(shí)間步長流入和流出該控制體的液態(tài)水質(zhì)量;mimp為當(dāng)前時(shí)間步長撞擊到控制體表面的液態(tài)水質(zhì)量;mevap為當(dāng)前時(shí)間步長控制體內(nèi)液態(tài)水蒸發(fā)的質(zhì)量。

當(dāng)前時(shí)間步長控制體內(nèi)液態(tài)水膜質(zhì)量為

假設(shè)水滴均勻撞擊在控制體表面,當(dāng)前時(shí)間步長撞擊到控制體表面的液態(tài)水質(zhì)量為

式中:v∞為遠(yuǎn)場(chǎng)來流速度;LWC為液態(tài)水含量;β為局部水收集系數(shù)。

當(dāng)前時(shí)間步長流出該控制體的液態(tài)水質(zhì)量為

式中:ln為流出面的底邊長度。

計(jì)算當(dāng)前時(shí)間步長控制體內(nèi)蒸發(fā)水的質(zhì)量,其中,蒸發(fā)速率可通過Chilton-Colburn比擬理論確定,利用對(duì)流換熱系數(shù)hc確定相應(yīng)位置的對(duì)流傳質(zhì)系數(shù)kc,如下:

式中:Sc為Schimidt數(shù);hc為對(duì)流換熱系數(shù);kc為對(duì)流傳質(zhì)系數(shù);cp為定壓比熱;Pr為普朗特?cái)?shù),從而可推導(dǎo)蒸發(fā)速率為

式中:cp,air為空氣的定壓比熱;MWwater和MWair分別為水和空氣的分子量;pv,sat為飽和蒸氣壓;pT和TT分別為總壓和總溫;T為當(dāng)前時(shí)間步長控制體溫度;p為壁面控制體絕對(duì)壓力;rh為相對(duì)濕度;在水膜表面該值為1;pv,∞,sat為入口飽和蒸氣壓;p∞為入口環(huán)境壓力。

同時(shí)可定義凍結(jié)系數(shù)f,是指在微元表面上實(shí)際結(jié)冰量和控制體內(nèi)所有水之比:

1.3 能量守恒方程

圖3為控制體內(nèi)的能量守恒示意圖,設(shè)流入能量為正。其中,流入控制體的能量包括水滴撞擊帶來的能量,從鄰近控制體流入的水帶來的能量和控制體內(nèi)結(jié)冰釋放的能量,流出控制體的能量包括對(duì)流換熱、導(dǎo)熱、熱輻射帶走的能量,蒸發(fā)或升華帶走的能量和流出控制體的水帶走的能量[7]。由于導(dǎo)熱和熱輻射對(duì)能量傳遞的影響較小,一般可忽略。

根據(jù)圖3可得能量守恒方程為

圖3 控制體能量守恒示意圖Fig.3 Schematic diagram of energy conservation of control volume

式中:cp,w為水的定壓比熱容;T為當(dāng)前時(shí)間步長控制體的溫度;imp為單位時(shí)間撞擊水?dāng)y帶的能量;evap為單位時(shí)間蒸發(fā)水帶走的能量;ice為單位時(shí)間所結(jié)冰層攜帶的能量;c為單位時(shí)間對(duì)流換熱帶走的能量。

式(16)可表示為

參考質(zhì)量守恒方程,對(duì)式(17)進(jìn)行離散求解,可得

式中:Told為上一時(shí)間步長控制體的溫度;等號(hào)左邊第1項(xiàng)對(duì)應(yīng)的是控制體內(nèi)水膜自身攜帶的能量隨時(shí)間的變化;第2項(xiàng)表示在x1和x2方向控制體內(nèi)所有的流入流出的液態(tài)水?dāng)y帶的能量。從而,可將每個(gè)時(shí)間步長結(jié)冰計(jì)算中的能量守恒方程表示為

式中:Q和Q分別為上一時(shí)間步長計(jì)算結(jié)束后控制體內(nèi)存在的水膜和冰層自身攜帶的能量;Qw和Qice分別為當(dāng)前時(shí)間步長控制體內(nèi)存在的水膜和冰層自身攜帶的能量;Qimp為撞擊水?dāng)y帶的能量;Qin和Qout分別為流入和流出該控制體的水?dāng)y帶的能量;Qevap為水膜蒸發(fā)帶走能量;Qc為對(duì)流換熱帶走的能量。

上一時(shí)間步長計(jì)算結(jié)束后控制體內(nèi)存在的水膜和冰層自身攜帶的能量為

式中:cp,ice為冰的定壓比熱容;Lfus為結(jié)冰時(shí)的固化潛熱;T0=273.15K。

當(dāng)前時(shí)間步長控制體內(nèi)存在的水膜和冰層自身攜帶的能量為

撞擊水?dāng)y帶的能量為

式中:T∞為遠(yuǎn)場(chǎng)來流溫度。

流出該控制體的水?dāng)y帶的能量為

水膜蒸發(fā)帶走的能量Qevap,在不同結(jié)冰情況下采用不同的計(jì)算公式。

若沒有發(fā)生結(jié)冰,即控制體內(nèi)是液態(tài)水膜時(shí):

若發(fā)生部分結(jié)冰,即結(jié)明冰情況下,控制體內(nèi)溫度處于相變平衡溫度時(shí):

若全部結(jié)冰,即結(jié)霜冰情況下:

式中:Levap和Lsub為凝固點(diǎn)溫度時(shí)的蒸發(fā)潛熱和升華潛熱。

對(duì)流換熱帶走的能量:

令Q=Qin+Qimp-Qevap-Qout-Qc,聯(lián)立式(20)～式(30),對(duì)能量守恒方程進(jìn)行求解,可得到結(jié)冰冰層質(zhì)量為

2 求解方法

通過動(dòng)量方程、質(zhì)量和能量守恒方程及引入的凍結(jié)系數(shù),聯(lián)立方程式(6)、式(10)、式(15)、式(31)可以對(duì)結(jié)冰數(shù)學(xué)模型進(jìn)行求解。

基于Fluent計(jì)算得到的外部空氣-水滴流場(chǎng)數(shù)據(jù)[22],利用Fluent提供的用戶自定義函數(shù)(user defined function,UDF),求解旋轉(zhuǎn)結(jié)冰數(shù)學(xué)模型,發(fā)展了一種循環(huán)遍歷控制體的非穩(wěn)態(tài)迭代求解方法。具體步驟如下:

步驟1在上一時(shí)間步水膜與冰厚的基礎(chǔ)上(初始時(shí)刻都為0),對(duì)動(dòng)量方程式(6)迭代求解,得到所有微元控制體內(nèi)水膜的平均速度和相應(yīng)的流動(dòng)方向。

步驟2根據(jù)水膜的流動(dòng)速度和方向,循環(huán)全部控制體,找到邊界及駐點(diǎn)位置作為計(jì)算的開端。

步驟3從駐點(diǎn)處開始對(duì)控制體展開計(jì)算,假設(shè)控制體內(nèi)初始溫度為0℃,并全部結(jié)冰,計(jì)算控制體內(nèi)的質(zhì)量和能量守恒方程,從而得到凍結(jié)系數(shù)f,對(duì)凍結(jié)系數(shù)進(jìn)行判斷,得到其結(jié)冰類型(明冰、霜冰、不結(jié)冰),根據(jù)結(jié)冰類型更新凍結(jié)系數(shù)、溫度、結(jié)冰量、溢流水量、水膜厚度等。利用更新的結(jié)果迭代計(jì)算動(dòng)量方程、質(zhì)量和能量守恒方程,直至收斂,該控制體計(jì)算完成,將當(dāng)前控制體溢流水流出量作為輸入值,賦值給臨近控制體,鄰近控制體流入的溢流水質(zhì)量可知。

步驟4繼續(xù)計(jì)算鄰近控制體,直至全部控制體計(jì)算完成,并利用更新后的水膜厚度,對(duì)全部控制體內(nèi)水膜速度進(jìn)行更新,返回步驟1開始計(jì)算,直至完成收斂判斷,則當(dāng)前時(shí)間步長計(jì)算完成,輸出當(dāng)前時(shí)間步長冰形。

步驟5進(jìn)入下一時(shí)間步長結(jié)冰計(jì)算,直至所設(shè)定結(jié)冰時(shí)間全部計(jì)算完成。

上述步驟的具體流程如圖4所示。

圖4 旋轉(zhuǎn)表面結(jié)冰求解流程Fig.4 Flow chart for solving icing on rotating surfaces

其中步驟3和4為結(jié)冰計(jì)算關(guān)鍵,對(duì)其進(jìn)行進(jìn)一步詳細(xì)說明:

1)假設(shè)控制體初始結(jié)冰狀態(tài)。假設(shè)初始控制體溫度為0℃,控制體內(nèi)結(jié)明冰,且水全部發(fā)生凍結(jié),可知mout=0,mw=0,Qout=0,Qw=0,且m和m對(duì)應(yīng)上一時(shí)間步長的液態(tài)水量和結(jié)冰量,為已知量。

通過上述假設(shè)可以對(duì)式(31)求解,得到結(jié)冰量,從而可根據(jù)式(15)得到凍結(jié)系數(shù)f。

2)對(duì)凍結(jié)系數(shù)f進(jìn)行判斷。

①當(dāng)f≥1時(shí),說明此時(shí)控制體結(jié)霜冰,控制體內(nèi)水全部發(fā)生凍結(jié),沒有溢流水流出,控制體內(nèi)不存在液態(tài)水,令f=1,則可得

②當(dāng)0＜f＜1時(shí),說明此時(shí)控制體結(jié)明冰,部分水發(fā)生凍結(jié),控制體溫度處于相變平衡溫度。控制體內(nèi)結(jié)冰量由式(31)可知,且聯(lián)立式(10)、式(11)、式(13)可得

聯(lián)立式(33)和式(34)可得水膜厚度為

觀察式(6)、式(35)可知,水膜流動(dòng)速度和水膜厚度相互影響,相互耦合,無法單獨(dú)求解,需要聯(lián)立求解。通過Jacobi迭代法和Gauss-Seidel迭代法,在對(duì)單一控制體迭代計(jì)算時(shí),將該控制體上一迭代步驟計(jì)算結(jié)果中的水膜厚度和水膜速度帶入當(dāng)前迭代步驟,循環(huán)迭代求解式(6)得到水膜速度用于判定流入流出水量和方向,再通過式(35)求出當(dāng)前迭代步驟的水膜厚度,與水膜速度一起用于下一迭代過程,直至水膜流動(dòng)速度和水膜厚度達(dá)到收斂,當(dāng)前控制體計(jì)算完成,可用于結(jié)冰計(jì)算。

在全部控制體結(jié)冰計(jì)算完成后,利用更新后的水膜厚度和水膜速度,從步驟1開始重新判定流動(dòng)速度和流動(dòng)方向,再次對(duì)所有控制體進(jìn)行結(jié)冰計(jì)算,直至前后2次計(jì)算誤差達(dá)到一定程度,完成收斂。通過計(jì)算單一控制體內(nèi)收斂和全部控制體收斂,在個(gè)體-整體雙重收斂的作用下,完成對(duì)水膜流動(dòng)速度和水膜厚度的計(jì)算。

通過迭代計(jì)算得到水膜厚度和水膜流動(dòng)速度后,即可得明冰狀態(tài)下的結(jié)冰參數(shù):

③當(dāng)f=0時(shí),說明此時(shí)控制體內(nèi)不結(jié)冰,只存在液態(tài)水膜,即mice=0,且由式(35)可知水膜厚度為

聯(lián)立式(6)和式(37),采用明冰計(jì)算時(shí)的迭代求解方法,可得水膜厚度和流動(dòng)速度,從而得到液態(tài)水狀態(tài)下相應(yīng)的結(jié)冰參數(shù):

3)通過對(duì)f的判斷可以得到控制體內(nèi)的結(jié)冰狀態(tài),從而計(jì)算得到當(dāng)前時(shí)間步長控制體的溢流水流出量,水膜厚度和結(jié)冰量,然后將當(dāng)前控制體溢流水流出量作為輸入值,賦值給臨近控制體。

4)進(jìn)行下一控制體計(jì)算,直至全部控制體計(jì)算完成,并利用更新后的水膜厚度,對(duì)控制體內(nèi)水膜速度進(jìn)行更新,返回步驟1重新計(jì)算,直至完成收斂,則當(dāng)前時(shí)間步長計(jì)算完成。

3 仿真計(jì)算分析

針對(duì)飛機(jī)上的所有旋轉(zhuǎn)部件結(jié)冰進(jìn)行研究,分析旋轉(zhuǎn)部件的結(jié)冰特性,為便于仿真模型擴(kuò)展到所有旋轉(zhuǎn)部件,參考文獻(xiàn)[17]的處理方式,將簡化后的等截面圓柱作為結(jié)冰壁面,其直徑為60mm,長度為200mm,對(duì)稱安置于直徑600mm的轉(zhuǎn)軸上,來流空氣沿轉(zhuǎn)軸軸向方向,如圖5所示。

圖5 旋轉(zhuǎn)圓柱模型示意圖Fig.5 Schematic diagram of rotating cylinder model

將入口來流邊界設(shè)為速度進(jìn)口,出口邊界設(shè)為壓力出口,結(jié)冰表面和轉(zhuǎn)軸均設(shè)為旋轉(zhuǎn)壁面,計(jì)算條件為:壓力101325Pa,迎角0°,結(jié)冰溫度263.15K,來流速度60m/s,轉(zhuǎn)速30rad/s,液態(tài)水含量(LWC)0.8g/m3,水滴直徑20μm,結(jié)冰時(shí)間420s。取時(shí)間步長0.02s與0.05s進(jìn)行驗(yàn)證計(jì)算,對(duì)比結(jié)果發(fā)現(xiàn)最終結(jié)冰冰形和水膜相近,本文取時(shí)間步長0.05s進(jìn)行計(jì)算,所得圓柱三維結(jié)冰冰形如圖6所示,圓柱中間截面上的二維冰形和水膜厚度隨時(shí)間變化如圖7和圖8所示。

圖6 旋轉(zhuǎn)圓柱三維冰形Fig.6 Three-dimensional ice shape of rotating cylinder

圖7 結(jié)冰冰形隨時(shí)間變化Fig.7 Ice shape changing over time

圖8 水膜厚度隨時(shí)間變化Fig.8 Water film thickness changing over time

由圖7和圖8可知,隨著結(jié)冰時(shí)間的增長,不斷有撞擊水撞擊到結(jié)冰壁面上,結(jié)冰表面水收集量逐漸增多,導(dǎo)致圓柱表面結(jié)冰厚度逐漸增加,當(dāng)結(jié)冰環(huán)境溫度不變時(shí),結(jié)冰能力基本不變,只有厚度隨時(shí)間發(fā)生明顯變化;且隨著時(shí)間的增長,水膜厚度和覆蓋范圍變化很小,水膜達(dá)到一定結(jié)冰時(shí)間后基本不變。同時(shí)由于旋轉(zhuǎn)的作用,圓柱表面上的結(jié)冰區(qū)域和水膜覆蓋區(qū)域均向圓柱下表面(-y方向)發(fā)生偏移,下表面的結(jié)冰區(qū)域和水膜覆蓋范圍大于上表面。

水膜厚度雖然整體變化趨勢(shì)符合預(yù)期,但是在駐點(diǎn)處數(shù)值偏大,通過對(duì)數(shù)據(jù)進(jìn)行分析,推測(cè)由于駐點(diǎn)處速度趨近于0,水膜在駐點(diǎn)處基本不發(fā)生溢流,液態(tài)水在該點(diǎn)產(chǎn)生積聚,導(dǎo)致駐點(diǎn)處水膜厚度偏大?？稍谥蟮难芯恐屑尤雺毫Φ挠绊?當(dāng)剪切力趨于0時(shí),將壓差作為駐點(diǎn)處水膜流動(dòng)的主要驅(qū)動(dòng)力,對(duì)模型算法進(jìn)一步改進(jìn)。

由于旋轉(zhuǎn)表面結(jié)冰實(shí)驗(yàn)研究匱乏,缺乏相關(guān)實(shí)驗(yàn)數(shù)據(jù)用以驗(yàn)證數(shù)值模擬的準(zhǔn)確性,因此為了驗(yàn)證所提模型的可靠性,將該模型算法的結(jié)冰計(jì)算結(jié)果與FENSAP軟件獲得的結(jié)冰冰形進(jìn)行對(duì)比,來驗(yàn)證所提模型的合理性。在轉(zhuǎn)速為30rad/s,結(jié)冰時(shí)間420s的工況下開展計(jì)算,將最終冰形結(jié)果進(jìn)行對(duì)比,如圖9和圖10所示。

觀察圖9和圖10的對(duì)比可知,兩者的結(jié)冰厚度和結(jié)冰范圍基本相同,其中本文所搭建的結(jié)冰模型的結(jié)冰范圍略微偏小;兩者的冰形變化趨勢(shì)也基本一致,在旋轉(zhuǎn)作用下同樣向圓柱下表面發(fā)生了偏移;水膜厚度變化曲線基本吻合,只有在駐點(diǎn)處由于速度趨近于0,使得結(jié)冰模型計(jì)算精度不夠,水膜厚度數(shù)值偏大,其余點(diǎn)處吻合良好。本文模型與FENSAP計(jì)算結(jié)果對(duì)比分析,證明本文模型的結(jié)冰冰形和水膜溢流計(jì)算結(jié)果合理,所提結(jié)冰模型求解方法正確。

圖9 本文模型和FENSAP冰形結(jié)果對(duì)比Fig.9 Comparison of ice shape between the proposed model and FENSAP

圖10 本文模型和FENSAP水膜厚度對(duì)比Fig.10 Comparison of water film thickness between the proposed model and FENSAP

通過改變轉(zhuǎn)速,研究轉(zhuǎn)速對(duì)旋轉(zhuǎn)圓柱表面結(jié)冰和水膜流動(dòng)的影響,將轉(zhuǎn)速分別設(shè)為0,30,60rad/s,其二維中間截面的結(jié)冰冰形和水膜厚度如圖11和圖12所示。

根據(jù)圖11和圖12比較3種不同轉(zhuǎn)速結(jié)果可知,隨著轉(zhuǎn)速的增加,水膜受到的慣性力和空氣剪切力等逐漸增加,圓柱表面液態(tài)水的溢流能力逐漸增強(qiáng),結(jié)冰范圍和水膜覆蓋范圍向圓柱下表面(-y方向)出現(xiàn)明顯偏移,轉(zhuǎn)速越大,偏移現(xiàn)象愈加明顯。且隨著轉(zhuǎn)速的增加,冰層厚度基本不變,水膜厚度略有減小。

圖11 結(jié)冰冰形隨轉(zhuǎn)速變化Fig.11 Ice shape changing with speed

圖12 水膜厚度隨轉(zhuǎn)速變化Fig.12 Water film thickness changing with speed

通過改變撞擊水滴的直徑,研究其對(duì)旋轉(zhuǎn)圓柱表面結(jié)冰和水膜流動(dòng)的影響,撞擊水滴直徑分別設(shè)為20,25,30μm,其二維中間截面的結(jié)冰冰形和水膜厚度如圖13和圖14所示。

圖13 結(jié)冰冰形隨水滴直徑變化Fig.13 Ice shape changing with droplet diameter

圖14 水膜厚度隨水滴直徑變化Fig.14 Water film thickness changing with droplet diameter

通過圖13和圖14比較3種不同撞擊水滴直徑的計(jì)算結(jié)果可知,隨著水滴直徑的增加,水滴慣性更大,不易受到氣流的作用而繞過迎風(fēng)表面,更容易和圓柱表面發(fā)生碰撞,水滴撞擊范圍和撞擊量相應(yīng)增加,同時(shí)由于當(dāng)前環(huán)境結(jié)冰能力恒定,冰層厚度隨水滴直徑變化較小,未凝結(jié)的液態(tài)水相應(yīng)增多,會(huì)在空氣剪切力和離心力等的作用下發(fā)生溢流,因此結(jié)冰范圍和水膜覆蓋范圍隨水滴直徑增長而逐漸增加,水膜厚度也逐漸增大。

通過改變液態(tài)水含量,研究其對(duì)旋轉(zhuǎn)圓柱表面結(jié)冰和水膜流動(dòng)的影響,液態(tài)水含量分別設(shè)為0.6,0.8,1.0g/m3,其二維中間截面的結(jié)冰冰形和水膜厚度如圖15和圖16所示。

由圖15和圖16比較3種不同液態(tài)水含量的仿真計(jì)算結(jié)果可知,隨著液態(tài)水含量的增長,相同時(shí)間撞擊到迎風(fēng)表面的水滴增多,導(dǎo)致圓柱表面的液態(tài)水逐漸增多,同時(shí)由于當(dāng)前結(jié)冰環(huán)境的結(jié)冰能力恒定,結(jié)冰厚度只發(fā)生了小幅度增長,可溢流的未凝結(jié)液態(tài)水逐漸增多,水膜厚度也相應(yīng)增加,在空氣剪切力和離心力等的作用下,未凝結(jié)水膜發(fā)生溢流,其覆蓋范圍明顯變大,結(jié)冰范圍也明顯增加。

圖15 結(jié)冰冰形隨液態(tài)水含量變化Fig.15 Ice shape changing with liquid water content

圖16 水膜厚度隨液態(tài)水含量變化Fig.16 Water film thickness changing with liquid water content

4 結(jié) 論

考慮旋轉(zhuǎn)部件表面的水膜流動(dòng)受到空氣剪切力與慣性力影響,建立了旋轉(zhuǎn)表面的結(jié)冰熱力學(xué)模型,通過個(gè)體-整體雙重收斂的方式進(jìn)行非穩(wěn)態(tài)迭代求解。以旋轉(zhuǎn)圓柱表面為研究對(duì)象,分析各影響因素對(duì)結(jié)冰和水膜流動(dòng)的作用,得出以下結(jié)論。

1)所提模型的非穩(wěn)態(tài)求解方法經(jīng)過驗(yàn)證,可以有效對(duì)三維旋轉(zhuǎn)部件表面的結(jié)冰冰形和水膜流動(dòng)進(jìn)行模擬,在駐點(diǎn)處水膜計(jì)算尚存在不足,可在之后的研究中對(duì)模型算法進(jìn)一步改進(jìn)。

2)由于旋轉(zhuǎn)的作用,圓柱表面的結(jié)冰區(qū)域和水膜覆蓋范圍會(huì)發(fā)生偏移,且隨著結(jié)冰時(shí)間的增長,旋轉(zhuǎn)圓柱表面的結(jié)冰厚度也隨之增加,其結(jié)冰區(qū)域和水膜溢流受結(jié)冰時(shí)間影響較小。

3)隨著轉(zhuǎn)速的增加,旋轉(zhuǎn)圓柱表面結(jié)冰區(qū)域和水膜覆蓋范圍的偏移現(xiàn)象愈加明顯,冰層厚度基本不變,水膜厚度略微減小。

4)隨著水滴直徑的增加,旋轉(zhuǎn)圓柱表面的冰層厚度變化較小,其結(jié)冰區(qū)域和水膜覆蓋范圍隨水滴直徑增長逐漸增加,水膜厚度也逐漸增大。

5)隨著液態(tài)水含量的增長,旋轉(zhuǎn)圓柱表面的結(jié)冰厚度出現(xiàn)小幅增長,水膜厚度相應(yīng)增加,結(jié)冰范圍和水膜覆蓋范圍也明顯變大。