飛機除冰沖擊射流多相流耦合模型研究

龔 淼,申遠航

(中國民航大學 航空工程學院, 天津 300300)

0 引言

飛機除冰效果和效率是影響航班起飛安全和準點率的重要環節。近年來我國冬季冰雪天氣增多,降雪范圍增大,除冰作業壓力逐年增大。目前民航局已經在北京大興機場和首都機場推廣了慢車除冰模式,對提升除冰作業效率產生了積極的影響。在慢車除冰作業中,除冰沖擊射流在飛機表面形成的速度場與壓力場對除冰效果有很大影響,研究沖擊射流沖擊在飛機表面的壓強與流速大小可為縮短慢車除冰融冰時間和優化作業參數提供重要參考。由于飛機除冰研究具有尺寸大、耦合的物理場數目多、涉及研究范圍廣等難點,很多國內外學者都使用數值模擬法來研究除冰液沖擊射流。因此將除冰液沖擊射流的多相流模型進行對比分析,選擇出適合研究飛機除冰沖射流的多相流模型可以有效提高研究效率,在飛機除冰研究中具有重要的工程意義。

國內外學者用來研究沖擊射流的平臺包括CFX、Ansys、OpenFOAM、Comsol等。Castillo[1]使用CFX對大壩的溢流瀑布沖擊射流進行了研究,得到了大壩自由落體沖擊射流的空氣夾帶程度與射流厚度、破裂形狀、沖擊底部壓力的關系,仿真結果與實驗結果具有良好的一致性。此外還發現空氣夾帶程度對CFX的網格具有較高的敏感性,而改變仿真中的湍流模型只會對沖擊射流滯點出的速度壓力值有影響。Faiza等[2]使用了CFX模塊研究了湍流液體射流流入靜止水池中形成的流場結構,通過建立兩相流-氣泡流模型,使用粒子圖像測速法監測水池表面以下區域的時均速度場。此外還用相間動量傳遞模型預測了下落流體的平均流量,預測結果與實驗值具有較高的一致性。Bel等[3]使用能量方法預測氣泡爆炸的影響,如果使用微射流來估計氣泡爆炸產生的相關壓力場,得到的結果更容易被PIF(壓力強度函數)和EPF(侵蝕功率函數)響應。Mahdavi等[4]使用Ansys-fluent研究了納米射流噴射到圓形轉盤之后的傳熱特性與流動狀態,通過流體體積模型將兩相間張力的作用體現出來,并且發現納米流體的體積分數越大,傳熱所需的扭矩就越大。Wang等[5]利用OpenFOAM模擬了深海起重機控制有效系統載荷受水動力沖擊的場景,其中6 DOF剛體運動求解器約束條件的引入保證了結果的準確性,研究發現波粒速度和波斜率是圓柱體上壓力不對稱的基本因素。周洪利等[6]利用OpenFOAM研究了豎直管內復雜的氣固液三相流動問題,提出了有效的水合物顆粒聚并破碎模型組合,發現了管內壓力損失梯度與流速的變化關系。Islas等[7]使用OpenFOAM模擬了粉塵爆炸中粉塵云的變化情況。該模型預測了粉塵云中雙渦流模式的形成,使用較大的顆?？梢栽黾宇w粒慣性,改善粉塵云的均勻性。Elaissi等[8]通過Comsol建立了一個湍流等離子噴射模型用來幫助合成硅納米粉末。研究考慮了納米粒子通過對流、擴散和熱泳的運輸,給出了等離子射流流場以及粒子的軌跡和熱傳導過程。周忠會等[9]基于Comsol研究了海底沉積物的熱流場分布,通過近似原位的電導熱加熱實驗,給沉積物模型加熱并且探究了加熱后的沉積物熱流場分布,這項研究對預測和評估海底天然氣水合物資源量有重要意義。Comsol支持多種求解方法包括有限元法、邊界元法、VOF法等,可以有效計算多物理場耦合的復雜流動模型,比較適合用來研究飛機除冰過程中的沖擊射流流場參數變化。

在沖擊射流方面,大量研究結果表明圓形出口的沖擊射流滿足射流的一般特性。根據沖擊射流的流場形狀和速度可以將沖擊射流分為3個明顯的部分:自由射流區、沖擊區和壁面射流區。在自由射流區,射流擴展半徑隨射流中心軸線長度增長而增大。在射流沖擊區,射流與沖擊壁面接觸,發生明顯彎曲,存在很大的壓力梯度。在壁面射流區,流動速度方向發生顯著變化與沖擊壁面垂直[10]。貝爾陶斯和拉賈拉南設計了大量實驗,結合量綱分析發現了沖擊射流滯點壓力與出射速度、射流高度、射流半徑[11]的定量關系。Zhou等[12]研究了泵射流水動力性能、噪聲性能,涉及空化侵蝕和尖端間隙渦流特性,轉子與定子及尾流場的相互作用。研究發現,緊湊的結構和復雜的內外部流場使泵獲得了更高的性能。

本文以飛機除冰噴灑作業為研究對象,基于二維的N-S方程,根據Comsol中的多相流模型模擬了飛機除冰作業中的除冰沖擊射流,得到了沖擊射流在沖擊壁面上的速度場與壓力場分布,并且根據速度與壓力參數確定了混合物模型和相傳遞混合模型更適合做除冰沖擊射流的仿真研究,為飛機除冰控制系統的研究工作提供了參考。

1 沖擊射流幾何及數值模擬

1.1 射流沖擊幾何模型

本研究以飛機除冰噴灑作業為研究對象,如圖1(a)所示。為了更加明顯地對比出不同多相流模型對模擬飛機除冰沖擊射流的模擬結果,將實際沖擊射流模型簡化,選取圓柱形噴槍的中心截面為研究平面,模擬垂直沖擊角度下所取二維研究平面內氣液兩相流場的速度與壓力分布情況,如圖1(b)所示。在Comsol中建立沖擊射流平面模型,根據除冰車噴嘴的實際參數確定模型尺寸,其中窄邊為除冰槍出口,寬邊為除冰槍內壁。為匹配實際除冰過程,選取1 m×1 m的沖擊空間,計算域如圖1(c)所示。根據行業標準AHM975規定,飛機除冰車噴槍出口溫度范圍為60～85 ℃,除冰車噴灑壓力約為9 bar,除冰作業時噴灑流量約為200 L/min,因此設置噴嘴直徑0.012 m,長0.05 m,數值模擬的流場入口為ab,出口選定ce和df,沖擊壁面為ef,流體屬性設置為標準大氣壓下85 ℃水的物性參數。

圖1 飛機除冰車噴嘴及計算域劃分圖

1.2 多相流模型的控制方程

本研究以相場模型、相傳遞混合模型、混合物模型、水平集模型、氣泡流模型這5個多相流模型為基礎,分別與Realizablek-ε[13]湍流模型耦合,添加穩態研究,探究穩態情形下的沖擊射流速度場與壓力場分布情況。

Realizablek-ε是標準k-ε[14]模型經過修正之后得出的一種湍流模型,更加適合用來解決流速變化大,流場形狀復雜的流動場景。Realizablek-ε的求解基于動量守恒的RANS方程、質量守恒的連續性方程、壁附近的壁函數以及k-ε輸運方程。

基于動量守恒的RANS方程:

ρ(u·▽)u= ▽·[ -ρI+K]+F

(1)

連續性方程:

▽(ρu)=0

(2)

壁面無滑移條件:

u·n= 0

(3)

式(1)中ρ是流體密度,u是速度場矢量,▽是二維平面下的拉普拉斯算子,F是包含重力在內的體積力。其中K是粘性應力張量。

K= (μ+μT)(▽u+ (▽u)T)-

(4)

湍流動能k的輸運方程:

(5)

湍流動能耗散率ε的輸運方程:

(6)

(7)

(8)

在式(4)和式(5)中ep是湍流耗散率;式(7)中μT是湍流粘度;Cμ是隨應變率變化的量,簡化計算可以取0.09;式(8)中Pk是生成項,Cε1取1.44,Cε2取1.92,σε取1.3,σk取1.0。Realizablek-ε與標準k-ε使用相同的默認初始條件和縮放比例。

相場模型使用式(9)追蹤相界面移動情況,其中λ是混合能量密度,γ移動參數。

(9)

相傳遞混合模型使用式(10)控制粘度,其中μd和μc是連續相和分散相的粘度。

(10)

水平集添加了式(11)跟蹤相界面,方程表征的是相界面的移動情況。其中γ是重新初始化參數,一般取1,ε是界面厚度控制參數,φ是水平集變量根據相的不同取1或者-1。

(11)

氣泡流模型結合了氣液兩相的動量方程和連續性方程,其動量方程為:

(12)

1.3 網格劃分及邊界條件

為了最大程度降低網格對不同流動模型仿真的影響,本研究的5個多相流模型都使用同等級的物理場控制的網格。改變網格最大單元尺寸,對比網格質量如表1所示。

表1 不同最大網格單元大小下的網格參數

由表1可以看出,改變網格尺寸,網格質量變化不大,保持在0.82～0.83,能較好地滿足計算精度。最終選取中間值,即總網格數為13 855,平均網格質量為0.828 1的網格劃分方案。網格劃分結果如圖2所示,網格質量分布如圖3所示。

圖2 計算域網格劃分圖

在圖3中,顏色越趨于藍色的網格質量越好,可以看出本研究的模型網格具有相對較好的網格質量。網格劃分部分使用了最常用的自由三角形網格,能夠在滿足計算收斂性的前提下減少計算量。為了防止邊界處的湍流變量和流體流動變量的不收斂,在湍流模型的壁面邊界進行了網格加密。此外,本研究旨在探究不同多相流模型在計算除冰液射流問題中的適用性,對不同模型采用了相同的網格劃分辦法,因此不再進行網格無關性驗證。

為保證計算準確性和穩定性設置了以下邊界條件:使用速度條件控制入口,根據行業標準AHM975將入口速度設置為30 m/s,法向流入;左右出口使用一個標準大氣壓力約束;液氣兩相流不可壓縮且包含重力作用;在噴嘴內壁和其他無熱流交換的壁面設置為無滑移條件;使用了P1+P1和一階線性的離散化方案,即用一階函數的辦法求解速度場和壓力場。此外還根據雷諾數指定了統一的中等湍流強度[15]。

2 結果對比分析

2.1 沖擊射流速度場結果對比分析

不同多相流模型的沖擊射流速度分布云圖如圖4所示。觀察速度云圖可以發現噴嘴入口處的速度最大為設置值30 m/s,從噴口處往下,射流在渦旋卷積作用下逐漸擴張,截面逐漸增大。速度沿中心流線逐漸減小,在射流到達沖擊壁面時速度急劇減小到0,而后在沖擊壁面上形成穩定的壁面射流[16]。

圖4 沖擊射流多相流模型速度分布云圖

然而觀察速度云圖可以發現相場模型和水平集模型在沖擊壁面處發生了偏轉,這與實際情況不符。研究2個模型的控制方程可知,相場模型和水平集模型都是基于場的模型,液相與氣相之間互不相容,在相鄰計算網格之間,相場函數φ[17]會在-1和1之間跳躍,故2種模型對相界面的描述較精確,但是難以準確模擬存在大量流體顆?；蛘呶F的分散型流場。綜上在圖4(a)和圖4(d)中射流發生偏轉是由于射流在破碎時形成了大量流體微團,上述2個模型無法準確描述這種分散相,導致射流兩側流場分布不對稱,隨著射流的不斷發展形成了偏轉的現象。

接下來探究沖擊壁面上的速度場分布情況,以沖擊壁面表面為研究對象,獲取不同多相流模型表面的速度值繪制曲線如圖5所示。

圖5 沖擊壁面速度大小分布曲線

不同模型形成的壁面沖擊射流具有大致相同的分布模式,在滯點[18]位置速度最小,相場和水平集模型由于其偏轉的程度最高,射流沖擊到壁面時與壁面存在夾角,動能沒有完全轉換為壓力能,最小速度分別為為0.36 m/s和0.46 m/s?；旌衔锬Ｐ?、氣泡流模型以及相傳遞混合模型在滯點處的動能都完全轉化為壓力能,最小速度為0。

壁面射流從滯點位置向兩側發展,壓力能逐漸轉換為動能,流動速度不斷上升。在出口處采用了標準大氣壓力約束,出口處是自由邊界,在壁面流動阻力的作用下,射流速度逐漸降低。不同模型在偏離滯點0.15～0.2 m的位置逐漸獲得最大流動速度。混合物模型的壁面最大流速為6.84 m/s,氣泡流模型的壁面最大流速為9.05 m/s,水平集模型的壁面最大流速為8.28 m/s,相場模型的壁面最大流速為6.84 m/s,相傳遞混合模型的壁面最大流速為6.67 m/s。其中氣泡流的最高速度比其他模型的最高速度都大,原因是氣泡流模型忽略了分散相的動量作用,流動過程中動能耗散更少,因此壁面射流的最高速度也最大。相場模型與水平集模型由于形成了偏轉,在滯點撞擊壁面時損失的能量更少,所以在壁面射流區獲得的最高速度大于相傳遞混合模型與混合物模型。

除冰作業的實際效果與滯點和最大速度點的位置有密切的關系。在今后的研究中,還將繼續探究慢車除冰作業的運動流場中,滯點位置的運動方式對除冰效果的影響。

2.2 沖擊射流壓力場結果對比分析

不同多相流模型的沖擊射流壓強分布圖如圖6所示。觀察壓強等值線圖,發現不同流動模型模擬出的壓強分布大致相同。在入口處流體靠壓力推動向前噴射,壓強達到最大值。隨著沖擊射流離開噴嘴,射流與大氣接觸,壓力能迅速釋放,流體壓強降低到大氣壓力。在滯點處,動能又迅速轉換為壓力能,滯點處獲得沖擊壁面的最大壓強。

圖6 沖擊射流多相流模型壓力分布圖

氣泡流的最大壓強為5.54×104Pa,明顯高于其他4個多相流模型的最大壓強。在2.1節中分析得出氣泡流模型在模擬沖擊射流的過程中忽略了分散相的動量作用,降低了能量損耗,導致其壁面射流最高速度遠遠高出其他4個模型,所以氣泡流模型在動能與壓力相互轉化的過程中也能獲得最大的壓強。

在沖擊射流發展段的兩側都形成了“眼狀”的等值線,從中心處向外逐漸減小到大氣壓。但是混合物和水平集模型的壓強等值線比其余3個模型稀疏,推測可能是其余3個模型模擬的射流左右兩側形成了渦旋,渦旋改變了周圍的壓強分布梯度。

通過檢測射流兩側的渦流強度來驗證猜想。取模型空間坐標系坐標為(-0.25 m,-0.5 m)處為監測點,獲取該點處不同多相流模型的渦流值。相場模型在該點處的渦流大小為37.742 s-1,相傳遞混合模型在該點處的渦流大小為6.86 s-1,氣泡流模型在該點處的渦流大小為20.521 s-1,水平集和混合物模型在該點未檢測到渦流。因此可以認為,在入口速度為30 m/s的情況下,相場、氣泡流、相傳遞混合模型在射流兩側形成了渦旋,造成了更大的壓強梯度。

2.3 沖擊射流的紊流系數大小對比分析

沖擊射流的紊流系數是射流結構的重要表征參數。一般來說出口處的射流紊流系數與紊流強度相關,紊流強度越大,紊流系數越大;一些關于垂直純水射流的研究發現當紊流系數取0.067時,實驗結果與理論計算結果具有較高的一致性[19]。為了進一步對比不同多相流模型對飛機除冰射流的模擬效果,將不同模型研究下的沖擊射流紊流強度進行對比。

在自由射流區,射流擴展半徑隨軸線長度呈線性增長,射流邊界與中心軸線形成張角α。一般定義射流張角α與紊流系數a的關系為:

tanα=ka

(13)

k是實驗系數,由實驗過程中的射流噴口形狀確定,本研究中由于出口相比于計算域而言較為狹窄,在工程中一般取2.44。因a是在理想情況下自由沖擊射流的紊流系數,本研究考慮了重力作用,故取距離噴口0.2 m以內的射流部分讀取其張角α。不同模型的射流邊界圖如圖7所示。

圖7 不同模型下的沖擊射流邊界圖

由圖7可以讀取不同模型下的沖擊射流邊界夾角,由大到小分別為8.5°、7.7°、7.25°、6.5°、6.25°。再根據式(13)得到紊流系數a,不同模型下的沖擊射流紊流系數如表2所示。

從表2中可以看出,混合物模型的紊流系數與實驗的紊流系數較為接近,水平集模型、相傳遞混合模型的計算結果與實驗值有一定差距,而相場模型、氣泡流模型的計算結果與實驗值相差最多,效果最差。

3 結論

1) 沖擊射流的截面沿著中心軸線呈錐狀發展,在滯點處流速衰減為0。不同模型在偏離滯點0.15～0.2 m的位置逐漸獲得最大流動速度。混合物、氣泡流、水平集、相場以及相傳遞混合模型在30 m/s的入射速度下于壁面獲得的最大流速分別為6.84、9.05、8.28、6.84、6.67 m/s。壁面的最大流速點到出口之間可以形成比較穩定的壁面射流。噴嘴入口處流體壓強最大,隨著沖擊射流離開噴嘴,流體壓強迅速降低到大氣壓力。在滯點處機械能相互轉化獲得沖擊壁面的最大壓強。

2) 相場模型和水平集模型在沖擊壁面處發生了偏轉,是由于沖擊射流在破碎時形成了大量流體微團,2個模型的相場函數φ會在-1和1之間突變,模型不能準確描述發生劇烈拓撲變化的流動,導致射流兩側流場分布不對稱,隨著射流的不斷發展出現了偏轉。

3) 通過沖擊射流的紊流系數大小對比可以發現,混合物模型、水平集模型、相傳遞混合模型的紊流系數與實驗值更加接近,模擬效果更好。綜合考慮模型模擬出的速度、壓力分布情況以及射流結構發現:對于飛機除冰射流以及類似的大流速、多流體微團類的流體分析,采用混合物模型和相傳遞混合模型可達到更加精確的研究結果。