郭 龍, 劉森云, 王 橋, 趙獻禮, 易 賢

(1.西北工業(yè)大學 航空學院, 西安 710072; 2.中國空氣動力研究與發(fā)展中心 結冰與防除冰重點實驗室, 四川 綿陽 621000)

0 引 言

飛機在穿越云層時可能遭遇過冷水滴，進而在迎風部位發(fā)生結冰現(xiàn)象。飛機結冰會改變其氣動外形，嚴重影響飛機的飛行性能和操作性能，威脅飛行安全[1]。1994年美國鷹航ATR72飛機發(fā)生空難，事故調(diào)查表明飛機遭遇了尺寸超過50 μm的過冷大水滴，發(fā)生了嚴重結冰，導致飛機操控失效而墜毀[2]。隨后，在1996-1997年冬季，NASA-Lewis研究中心聯(lián)合美國國家航空航天局(NASA)、美國國家大氣研究中心(NCAR)和美國聯(lián)邦航空管理局(FAA)進行了結冰云中水滴直徑測量的飛行試驗，在29次飛行試驗中出現(xiàn)了3次直徑超過50 μm的過冷大水滴[3]。加拿大國家研究委員會(NRC)開展了同樣的試驗，過冷大水滴的出現(xiàn)概率達到8%，水滴直徑覆蓋50～3000 μm[4]。此后，過冷大水滴環(huán)境下的結冰研究引起了人們的廣泛關注[5-7]。

結冰數(shù)值模擬中，有2個最基本的假設：一是水滴的剛性球假設，即水滴在運動過程中始終保持球形，不發(fā)生變形和破碎現(xiàn)象，水滴阻力采用經(jīng)典的球形阻力模型計算；二是水滴在碰撞固壁后不發(fā)生飛濺現(xiàn)象，即水滴撞擊量與結冰量守恒[8]。在常規(guī)飛機的飛行中，飛機周圍流場的壓力梯度一般不會導致小于50 μm的水滴產(chǎn)生變形和破碎等動力學特性[9]，因此針對小水滴,以上的假設是合理的。Gunn等[10]觀測了大水滴的降落過程，發(fā)現(xiàn)水滴由圓球形變成橢球形，并嚴重地影響水滴下降速度，這說明大水滴會發(fā)生變形現(xiàn)象，而變形會導致水滴阻力特性的變化。Tan等[11]通過實驗手段觀測到大水滴在靠近機翼的過程中會發(fā)生變形和破碎現(xiàn)象，在碰撞機翼時會發(fā)生飛濺現(xiàn)象。因此，前文介紹的2個基本假設不再適用于大水滴的結冰數(shù)值模擬。大水滴動力學特性的存在使得其結冰數(shù)值模擬更加復雜。

以往對水滴變形形態(tài)和阻力特性的研究通常采用水滴受重力作用而自由下落的實驗方式[12-14]，但是該實驗獲得的水滴變形量有限，不能充分展示水滴在飛機周圍流場中的變形。本文通過搭建水滴動力學特性實驗平臺，研制水滴發(fā)生器，利用高速相機記錄水滴在氣動力作用下的運動和變形過程，拓展水滴變形及阻力特性研究范圍，希望可為飛機大水滴結冰數(shù)值計算模型的修正提供實驗數(shù)據(jù)支持。

1 實驗設備和方法

1.1 實驗設備

整個實驗平臺在某0.3 m×0.2 m結冰風洞駐室的基礎上進行改造搭建，如圖1所示。實驗風道采用回流式設計，通過2個拐角將駐室、豎直風道、風機相連，形成回路。結冰風洞駐室配有制冷系統(tǒng)，能夠調(diào)節(jié)風道內(nèi)的氣流溫度；豎直風道和拐角段以透光率超過90%的有機玻璃制成。風道收縮段采用二元維多辛斯基曲線設計，收縮比為4；收縮段長0.2 m，其入口尺寸為0.3 m×0.3 m，出口尺寸為0.3 m×0.075 m。水滴的變形過程主要發(fā)生在收縮段，為避免曲面給拍攝結果帶來誤差，將高速相機安裝于二元平面?zhèn)?，采用鹵素燈背光方式照明。

圖1 實驗系統(tǒng)圖

根據(jù)Plateau-Rayleigh不穩(wěn)定的原理，水柱會自然地趨向不穩(wěn)定、發(fā)生振蕩，形成不同大小的水滴。壓電式等徑離散水滴發(fā)生器(見圖2)通過壓電陶瓷振動向水柱引入一個合理的初始擾動，提供水柱表面波形成的初始波長，該初始波不斷發(fā)展成為水柱振蕩的主波，當波的振幅接近水柱直徑時，波幅不能繼續(xù)增大，水柱將被打斷形成離散的等徑水滴，如圖3所示。水滴的大小可以通過水壓、孔口直徑、壓電片振動頻率和振幅聯(lián)合控制，實驗調(diào)試獲得的水滴直徑范圍為240～910 μm。較大直徑的水滴可以通過點滴器獲取，其水滴直徑取決于點滴針頭內(nèi)徑與剪切氣流的速度。

圖2 壓電式等徑離散水滴發(fā)生器

圖3 等徑離散水滴

1.2 實驗方法

收縮段內(nèi)氣流速度的分布是利用空速管對其入口和出口風速進行測量、再利用數(shù)值模擬的方法計算獲得。水滴的運動過程利用HG-100K高速相機進行拍攝，相機參數(shù)選取1024 pixel×1024 pixel，幀數(shù)為1500幀/s。拍攝前對水滴進行對焦，然后將標尺固定在水滴焦平面進行標尺拍攝，利用標尺建立單位像素與單位長度的關系。水滴圖片經(jīng)過處理后，利用圖片中的像素信息可以計算出水滴的各項所需參數(shù)。

利用MATLAB程序?qū)λ螆D片進行處理。處理過程大致如下：首先對圖片進行灰度處理，然后通過中值濾波和高斯濾波對圖片進行降噪與光順，利用Sobel算子進行水滴邊緣的初步檢測，再進行非極大抑制并用雙閾值算法最終完成水滴的邊緣檢測和連接，最后去除水滴內(nèi)部點進行重心位置檢測。各步驟處理結果如圖4所示。

圖4 水滴圖像處理

2 實驗結果和分析

2.1 收縮段速度分布

利用空速管測量不同電機轉(zhuǎn)速下收縮段入口與出口的氣流速度，測試結果見表1。以入口氣流速度作為邊界條件，出口氣流速度作為壓力邊界條件，利用Fluent軟件計算整個收縮段與實驗段內(nèi)的氣流速度分布情況。不同條件下風道收縮段與實驗段中心線氣流速度分布如圖5所示。S為收縮段入口到實驗段中心線的距離。將收縮段出口氣流速度的計算結果與實驗測試結果進行對比(見圖6)，計算結果與實驗結果吻合較好，驗證了計算方法的正確性和速度參數(shù)的可靠性。

表1 不同轉(zhuǎn)速下收縮段入口與出口氣流速度

圖5 收縮段和實驗段中心線氣流速度

圖6 收縮段出口氣流速度計算值與實測值對比圖

2.2 水滴的變形

由于表面張力的作用，水滴會自發(fā)地趨向于球形，以達到最小的表面積和表面能。對于穩(wěn)定的球形水滴，其對稱面受向內(nèi)的表面張力和由內(nèi)外壓差形成的向外的壓力作用，兩個力大小相等，方向相反，如圖7所示：

2πrσ=Δpπr2

(1)

即有：

(2)

式中,r為水滴半徑，Δp為水滴內(nèi)外壓差，σ為水滴表面張力系數(shù)，取值0.073 N/m。

圖7 穩(wěn)定球形水滴受力

由式(2)可知，水滴受到均勻穩(wěn)定的外部壓力作用時，能夠保持球形；水滴受到不均勻表面壓力作用時，表面受力將失衡，發(fā)生變形甚至破碎現(xiàn)象。水滴在氣流作用下加速運動時便是這種情況。在氣流作用下加速運動的水滴,其形狀主要由以下3個因素共同決定[9]：

(1) 由表面張力產(chǎn)生的附加壓強

ps=4σ/d

(3)

(2)水滴外部動壓

p外=ρaU2/2

(4)

(3)由加速度產(chǎn)生的內(nèi)部靜壓

p內(nèi)=ρwad

(5)

式中,d為水滴等效直徑，ρa為空氣的密度，U為氣流與水滴的速度差，ρw為水的密度，a為水滴的加速度。

無量綱數(shù)We表征慣性力與表面張力之比，其值為水滴外部動壓與附加壓強之比的8倍：

We=ρaU2d/σ=p外/8ps

(6)

無量綱數(shù)Bo表征重力與表面張力之比，其值為水滴內(nèi)部靜壓與附加壓強之比的4倍：

Bo=ρwad2/σ=p內(nèi)/4ps

(7)

可見，水滴的形狀與We、Bo有很大關系。

圖8為水滴在氣動力作用下變形的演化過程。水滴直徑為2.18 mm，收縮段入口和出口氣流速度分別為8.84 m/s和35.90 m/s。水滴由球形逐漸變扁，最終形成圓盤形，整個變化過程所用時間約為4.69 ms。

圖8 水滴變形演化過程

通過實驗觀察，水滴在變形過程中隨著We和Bo的增大，依次歷經(jīng)了4種典型的形狀(如圖9所示)：①圓球形，水滴各個方向曲率半徑相當；②橢球形，水滴迎風面和背風面曲率半徑增大，并大于側面曲率半徑；③半球形，水滴迎風面變?yōu)槠矫?，背風面仍為曲面，曲率半徑繼續(xù)增大；④圓盤形，水滴迎風面與背風面均呈平面。

利用水滴縱橫比E來表征水滴變形程度：

圖9 水滴變形形狀

(8)

式中，dh為水滴縱向厚度，dv為水滴橫向?qū)挾取?/p>

水滴縱橫比E與We的關系如圖10所示。隨著We增大，E呈線性減小。實驗中氣流和水滴經(jīng)過收縮段時都作加速運動，而水滴加速運動的動力主要來源于氣流，因此,氣流與水滴速度差會沿著收縮段增大，水滴外部動壓增大，We增大，水滴變形量增大，E減小。

圖10 水滴縱橫比與We的關系

水滴縱橫比E與Bo的關系如圖11所示。隨著Bo的增大，E減小，呈雙曲線變化。水滴收縮段加速運動時，加速度不斷增大，Bo增大，水滴變形量增大，E減小。

圖11 水滴縱橫比與Bo的關系

由式(6)和(7)可知，We只體現(xiàn)了水滴外部動壓與附加壓強對水滴變形的影響，Bo只體現(xiàn)了水滴內(nèi)部靜壓與附加壓強對水滴變形的影響，兩者均未完全考慮影響水滴形狀的3個因素。在此，結合We與Bo，引入新的無量綱參數(shù)Sn：

(9)

根據(jù)實驗數(shù)據(jù)，建立水滴縱橫比E與Sn的關聯(lián)式：

E=a1×e(-Sn/t1)+b1×e(-Sn/t2)+c1

(10)

如圖12所示，利用實驗數(shù)據(jù)，通過非線性擬合，得出a1、b1、c1、t1、t2的值分別為0.508 91、0.176 82、0.291 56、160.462 18、9.161 61。將其帶入式(8)可得E與Sn的關系式：

(11)

圖12 水滴縱橫比與Sn的關系

圖13為式(11)計算所得E與實驗測得E的對比圖。計算所得E與實驗測量值的偏差絕對值大部分都在10%以內(nèi)，說明Sn與E所建立的關系式能夠很好地預測水滴變形后的縱橫比。

圖13 水滴縱橫比的實驗測量值與計算值對比

2.3 變形對水滴阻力特性影響

阻力是物體在流體中作相對運動時產(chǎn)生的與相對運動方向相反的力。本實驗中水滴受氣流阻力和重力作用加速向下運動，水滴阻力系數(shù)定義為水滴所受到的阻力與氣流動壓和迎風面積之比。對水滴進行受力分析，利用牛頓第二定律可得：

(12)

式中，g為重力加速度，CD為水滴阻力系數(shù)，Dx為水滴迎風面直徑。

由式(12)可得出考慮水滴變形效應的CD計算表達式：

(13)

式中，ρd為水滴的密度。水滴在氣動力作用下會發(fā)生變形現(xiàn)象，水滴變形會增大水滴迎風面的面積，同時迎風面曲率半徑增大會促進水滴后緣更早地流動分離，這將導致水滴的阻力增大。引入考慮變形效應的水滴雷諾數(shù)Re，計算式如下：

(14)

式中，μa為空氣動力黏度。

根據(jù)實驗結果，可得CD隨Re變化的曲線，如圖14所示。由圖可知，在本實驗Re范圍內(nèi)，CD隨著Re的增大先略有減小后逐漸增大。根據(jù)實驗數(shù)據(jù)建立的CD擬合關系式如下：

CD=-170.48+206.65x-92.29x2+18x3-1.29x4

(15)

式中，x=log10Re。

圖14 水滴阻力系數(shù)隨Re變化

水滴變形主要歷經(jīng)了4種典型的形狀：圓球形、橢球形、半球形和圓盤形。圖15為圓球形阻力系數(shù)[15]、圓盤形阻力系數(shù)[16]和水滴阻力系數(shù)隨Re變化的對比圖。由圖可知，本文測得的CD與Gunn等測量自由下落水滴獲得的阻力系數(shù)非常接近，并且本文中水滴變形直至破碎，測得的CD很好地擴展到接近圓盤形阻力系數(shù)；水滴變形后阻力系數(shù)明顯增大，CD在Re≤500時與圓球形阻力系數(shù)吻合較好，當Re>500時CD開始偏離圓球形阻力系數(shù)，并不斷增大至圓盤形阻力系數(shù)，這與水滴變形形狀相對應。結合Clift等給出的球形阻力系數(shù)模型和本文測量擬合的水滴變形后的阻力系數(shù)模型，給出考慮變形的水滴阻力計算模型如下：

(16)

圖15 水滴阻力系數(shù)與圓球形、圓盤形阻力系數(shù)對比

3 結 論

采用實驗的手段對水滴加速運動過程中的變形現(xiàn)象及其對水滴阻力特性的影響進行了研究，結果表明：

(1) 水滴在加速過程中會發(fā)生變形現(xiàn)象，變形主要歷經(jīng)圓球形、橢球形、半球形和圓盤形4個過程。

(2) 隨著We和Bo的增大，水滴變形量增大，水滴縱橫比E分別呈線性和雙曲線性減小。

(3) 引入的無量綱數(shù)Sn與E的關聯(lián)式能夠很好地預測水滴變形量，實驗測量值與計算值偏差的絕對值大部分在10%以內(nèi)。

(4) 水滴變形會導致水滴阻力系數(shù)CD增大，其增大過程與其變形過程一致，由圓球形阻力系數(shù)向圓盤形阻力系數(shù)過渡。