過(guò)冷大水滴變形與破碎的影響因素

李維浩，易賢，李偉斌，王應(yīng)宇

中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000

飛機(jī)結(jié)冰是威脅飛行安全的主要因素之一[1-2]。長(zhǎng)期以來(lái)，研究者一直認(rèn)為飛機(jī)結(jié)冰是由平均容積直徑小于50 μm的過(guò)冷小水滴所致，因此相關(guān)結(jié)冰研究也主要針對(duì)過(guò)冷小水滴開展[3]。1994年ATR-72空難之后[4]，過(guò)冷大水滴(Supercooled Large Droplet, SLD)結(jié)冰問(wèn)題開始引起關(guān)注，關(guān)于SLD結(jié)冰的研究熱潮一直持續(xù)至今。研究表明，流場(chǎng)中的SLD在與物面相向運(yùn)動(dòng)的過(guò)程中，存在明顯的變形，并在一定條件下會(huì)發(fā)生破碎，這種水滴的變形與破碎現(xiàn)象會(huì)對(duì)結(jié)冰產(chǎn)生直接影響。因此，在進(jìn)行過(guò)冷小水滴結(jié)冰研究時(shí)關(guān)于水滴在運(yùn)動(dòng)過(guò)程中保持球體、不會(huì)變形和破碎的假設(shè)，并不適用于SLD結(jié)冰研究[5]。

針對(duì)SLD的變形與破碎問(wèn)題，國(guó)內(nèi)外做了大量的實(shí)驗(yàn)研究[6-8]，并提出了多個(gè)描述水滴變形與破碎的模型，這些實(shí)驗(yàn)都是在特定流場(chǎng)中進(jìn)行的，缺乏普適性，所以將這些經(jīng)驗(yàn)?zāi)Ｐ陀糜谀M水滴在機(jī)翼周圍復(fù)雜流場(chǎng)的變形與破碎缺乏準(zhǔn)確性[9]。針對(duì)這一不足，2005年Tan和Papadakis[10]提出用泰勒類比破碎(Taylor Analogy Breakup, TAB)模型[11]來(lái)描述水滴的變形與破碎過(guò)程，大大提高了SLD結(jié)冰數(shù)值計(jì)算的精度。雖然TAB模型只能模擬一種振蕩模式，而水滴破碎通常由多種振蕩模式共同作用導(dǎo)致，但目前還沒(méi)有可靠的破碎模型來(lái)替代[12]，因此TAB模型依然被廣泛應(yīng)用于SLD結(jié)冰計(jì)算中。目前，國(guó)內(nèi)外關(guān)于SLD變形與破碎數(shù)值計(jì)算方面的研究主要集中在變形與破碎對(duì)水滴收集系數(shù)和冰形的影響上[13-15]，缺少關(guān)于變形與破碎影響因素的系統(tǒng)研究。鑒于此，本文采用數(shù)值計(jì)算的手段，基于水滴運(yùn)動(dòng)計(jì)算的拉格朗日法和TAB模型，系統(tǒng)研究了水滴直徑、流場(chǎng)速度和繞流物體尺寸對(duì)水滴破碎與變形及其效應(yīng)的影響規(guī)律，在此基礎(chǔ)上，提出了一個(gè)表征水滴變形與破碎及相應(yīng)效應(yīng)的綜合影響參數(shù)。

1 數(shù)值計(jì)算方法

1.1 空氣流場(chǎng)

空氣流場(chǎng)通過(guò)求解低速黏性流場(chǎng)的時(shí)均Navier-Stokes方程獲得，采用SIMPLE系列算法求解[16]，控制方程的通用形式為

(1)

式中：t為時(shí)間；φ為輸運(yùn)變量；ρa(bǔ)為空氣密度；ua為空氣速度；Γφ為擴(kuò)散系數(shù)；qφ為源項(xiàng)。φ、Γφ和qφ取不同的值，可代表流場(chǎng)的連續(xù)性方程、動(dòng)量方程以及湍動(dòng)能、湍流耗散率等標(biāo)量的輸運(yùn)方程。式(1)中各項(xiàng)的物理意義、具體表達(dá)式和方程的求解方法，可詳見文獻(xiàn)[16]。

1.2 水滴軌跡

獲得空氣流場(chǎng)之后，本文采用拉格朗日法[17-18]進(jìn)行水滴軌跡計(jì)算。根據(jù)牛頓第二定律，水滴運(yùn)動(dòng)方程可以寫為

(2)

式中：xd為水滴位移；ρd為水滴密度；Vd為水滴體積；g為重力加速度；Ad為水滴的迎風(fēng)面積；CD為阻力系數(shù)；ua為當(dāng)?shù)亓鲌?chǎng)速度；ud為水滴速度。引入相對(duì)雷諾數(shù)：

(3)

式中：Deq為水滴當(dāng)量直徑；μa為空氣黏性系數(shù)。則式(2)可改寫為

(4)

以水滴的初始位置和流場(chǎng)的速度分布作為定解條件，水滴運(yùn)動(dòng)方程的求解屬于一階常微分方程的初值問(wèn)題，本文采用一階歐拉法進(jìn)行數(shù)值求解。

1.3 水滴變形與破碎

采用TAB模型[10-11]模擬水滴變形與破碎。TAB模型是假定水滴振蕩運(yùn)動(dòng)類似于彈簧-質(zhì)量-阻尼系統(tǒng)的受迫運(yùn)動(dòng)(如圖1所示)。

圖1 泰勒類比破碎模型Fig.1 Taylor analogy breakup model

描述彈簧-質(zhì)量-阻尼系統(tǒng)的受迫運(yùn)動(dòng)的數(shù)學(xué)方程可以寫為

(5)

式中：位移x是位于初始球體赤道處的點(diǎn)的位移；外力F是施加在水滴上的氣動(dòng)力；恢復(fù)系數(shù)k是由水滴表面張力引起的；阻尼系數(shù)f為水滴黏度。引入系數(shù)Cb、CF、Ck和Cd，且定義

Y=x/(Cbr)

(6)

(7)

(8)

(9)

式中：Y為水滴的無(wú)量綱變形量；r為水滴半徑。則式(5)可寫為

(10)

其中：

(11)

式(10)采用如下差分格式進(jìn)行離散求解：

(12)

式中：

(13)

當(dāng)位移超過(guò)初始球體半徑的一半，即Y>1時(shí)，認(rèn)為水滴破碎發(fā)生[10]，破碎后產(chǎn)生的子水滴直徑與母水滴直徑的關(guān)系表示為

(14)

其中：D為母水滴直徑；D32為子水滴粒徑分布的Sauter平均直徑[19]，定義為對(duì)應(yīng)尺寸水滴的體積與表面積之比等于所有水滴的總體積與總表面積之比；K為水滴的變形和振動(dòng)能量與基頻總能的比值，計(jì)算中K=10/3。子水滴直徑服從χ2分布[11]，即

(15)

式中：Γ(x)為伽馬函數(shù)；n為該分布的自由度，可取任意大于0的整數(shù)，子水滴直徑為自變量。由Sauter平均直徑的定義可求得n與D32的關(guān)系為

n=D32+2

(16)

子水滴的數(shù)量可由質(zhì)量守恒求得，且認(rèn)為子水滴是無(wú)變形、無(wú)振蕩的，并且具有與母水滴相同的速度。

1.4 水滴阻力

常見的飛機(jī)飛行速度范圍以及云層中經(jīng)常出現(xiàn)的水滴尺寸范圍，都超出了使用Stokes公式計(jì)算水滴黏性阻力的理論允許范圍，因此不能使用Stokes阻力計(jì)算公式，而使用低雷諾數(shù)下球體的嚴(yán)格阻力系數(shù)為[20]

(17)

水滴變形會(huì)引起阻力增加，阻力系數(shù)進(jìn)行以下修正[10]：

CD=CD,S(1+2.632Y)

(18)

1.5 數(shù)值計(jì)算流程

數(shù)值計(jì)算的流程可以總結(jié)為：① 輸入水滴初始的位置、直徑等信息，讀取之前已經(jīng)獲得的流場(chǎng)信息；② 判斷水滴所處網(wǎng)格位置，得到當(dāng)?shù)亓鲌?chǎng)速度，使用一階歐拉法迭代求解水滴軌跡，并且每迭代一步都判斷水滴是否撞擊到物體表面或者飛出流場(chǎng)；③ 使用TAB模型計(jì)算水滴的無(wú)量綱變形量Y，當(dāng)Y≥1時(shí)發(fā)生破碎，計(jì)算出破碎產(chǎn)生子水滴的平均直徑，采用相同方法計(jì)算子水滴的軌跡，并計(jì)算子水滴的粒徑分布和位置分布。

計(jì)算水滴軌跡及變形與破碎的具體流程如圖2 所示。

圖2 數(shù)值計(jì)算流程Fig.2 Workflow of numerical calculation

2 SLD變形與破碎特性計(jì)算分析

2.1 計(jì)算條件

SLD在機(jī)翼附近變形與破碎和水滴直徑以及水滴和空氣的相對(duì)速度有關(guān)，而相對(duì)速度主要決定于初始來(lái)流速度和物體的尺寸外形。對(duì)于NACA四位數(shù)字翼型，前緣半徑rL≈1.10c2/b，c為翼型厚度，b為弦長(zhǎng)。因此，本文選取弦長(zhǎng)均為1 m的NACA0008、NACA0012、NACA0015 3種厚度不同的對(duì)稱翼型，來(lái)流速度為90、120、150 m/s，水滴直徑為50、100、150、200、300 μm作為計(jì)算條件。

2.2 典型SLD流場(chǎng)區(qū)域劃分

圖3是直徑300 μm的水滴在速度150 m/s的流場(chǎng)下，在NACA0015翼型附近發(fā)生破碎前后的典型軌跡，圖3(b)是破碎點(diǎn)附近(圖3(a)的虛線圈)的放大圖。水滴破碎的原因是受到空氣剪切作用，當(dāng)水滴表面張力小于其受到的空氣剪切力時(shí)，發(fā)生破碎。當(dāng)水滴離翼型較遠(yuǎn)時(shí)，運(yùn)動(dòng)速度與流場(chǎng)一致，與空氣之間沒(méi)有相互剪切力，而當(dāng)水滴與翼型距離減小時(shí)，與空氣的速度差增大，在剪切力作用下，水滴先是變形，進(jìn)而破碎成若干子水滴繼續(xù)前進(jìn)，并碰撞在翼型的不同區(qū)域和位置。

從上面分析可知，包含SLD的典型流場(chǎng)可以劃分成以下幾個(gè)區(qū)域：① 水滴未撞擊區(qū)域，② 水滴撞擊但未破碎區(qū)域；③ 水滴待破碎區(qū)域；④ 水滴破碎區(qū)域。具體分布如圖4所示，圖中藍(lán)色區(qū)域的水滴不會(huì)撞擊在翼型表面，灰色區(qū)域水滴在撞擊翼型表面之前僅發(fā)生變形不會(huì)發(fā)生破碎，綠色區(qū)域水滴即將發(fā)生破碎，褐色區(qū)域是水滴破碎后存在子水滴的區(qū)域。在水滴待破碎區(qū)域，水滴變形量隨時(shí)間和空間推移而增加，越靠近翼型表面變形量越大，在兩區(qū)域交界處發(fā)生破碎。

圖3 水滴破碎前后的典型軌跡Fig.3 Typical trajectories of droplet before and after breakup

圖4 SLD結(jié)冰流場(chǎng)典型區(qū)域劃分Fig.4 Typical zoning of SLD icing flow field

2.3 SLD破碎區(qū)域面積

表1～表3分別給出了3種翼型流場(chǎng)在不同速度和水滴直徑條件下的破碎區(qū)域面積?？梢钥吹剑谒杏?jì)算狀態(tài)中，只有幾個(gè)速度高、水滴直徑大的狀態(tài)存在水滴破碎區(qū)域。水滴破碎區(qū)域面積受到流場(chǎng)速度、水滴直徑以及翼型前緣半徑的共同影響，速度越高、直徑越大、前緣半徑越大，水滴破碎區(qū)域面積越大。

表1 NACA0008翼型水滴破碎面積Table 1 Breakup area of NACA0008 airfoil mm2

表2 NACA0012翼型水滴破碎區(qū)域面積Table 2 Breakup area of NACA0012 airfoil mm2

表3 NACA0015翼型水滴破碎區(qū)域面積Table 3 Breakup area of NACA0015 airfoil mm2

2.4 水滴收集特性

圖5是有/無(wú)變形與破碎模型極限軌跡的對(duì)比(速度v=150 m/s，水滴直徑d=300 μm，NACA0015翼型)。圖中可以看出添加變形與破碎模型后撞擊極限變小，這是由于水滴變形引起阻力系數(shù)增大，水滴的跟隨性增強(qiáng)，軌跡的偏折程度會(huì)變大。圖6給出的是有/無(wú)變形與破碎模型的水滴收集系數(shù)β對(duì)比，橫坐標(biāo)是碰撞點(diǎn)離駐點(diǎn)的物面距離s與弦長(zhǎng)b的比值。其中圖6(a)顯示的是直徑50 μm的SLD在速度150 m/s的流場(chǎng)中NACA0015翼型上的水滴收集系數(shù)，雖然該流場(chǎng)沒(méi)有水滴破碎，但由于SLD在運(yùn)動(dòng)中會(huì)發(fā)生變形，變形導(dǎo)致水滴阻力系數(shù)增大，水滴運(yùn)動(dòng)軌跡與空氣流線一致性更高，從而導(dǎo)致局部水收集系數(shù)和撞擊極限都變小。圖6(b)顯示的是在相同流場(chǎng)下直徑300 μm的SLD的水滴收集系數(shù)，該流場(chǎng)包含有較大的水滴破碎區(qū)域，考慮變形與破碎之后，局部水收集系數(shù)和撞擊極限變化程度比圖6(a)狀態(tài)更大。

從圖6可以看出，不同的結(jié)冰條件，變形與破碎對(duì)水滴收集特性影響不一樣。圖7(a)給出的是對(duì)于NACA0012翼型，不同速度條件下變形與破碎對(duì)水滴總收集系數(shù)影響的曲線；圖7(b)顯示的是速度為120 m/s時(shí)，不同翼型條件下變形與破碎對(duì)水滴總收集系數(shù)影響的曲線，圖中的縱坐標(biāo)表示的是考慮變形與破碎之后水滴總收集系數(shù)變化的比例?？梢钥闯鏊沃睆健⒘鲌?chǎng)速度、機(jī)翼前緣半徑等與變形和破碎影響水滴總收集系數(shù)的程度呈正相關(guān)關(guān)系。速度為90 m/s時(shí)，影響程度較小，總收集系數(shù)變化不超過(guò)4%，速度增加至150 m/s之后，最大變化量可達(dá)到10%以上。

圖5 極限軌跡對(duì)比Fig.5 Comparison of limiting trajectories

圖6 局部水滴收集系數(shù)對(duì)比Fig.6 Comparison of local droplet collection coefficients

2.5 子水滴分布

圖8給出的是破碎后產(chǎn)生的子水滴的粒徑分布，母水滴是在同一流場(chǎng)狀態(tài)下(v=150 m/s,d=150 μm, NACA0012)，在遠(yuǎn)場(chǎng)的不同縱坐標(biāo)處釋放的，縱坐標(biāo)y分別為0和0.05?？梢钥闯觯纹扑楹?，會(huì)產(chǎn)生大量的子水滴，其分布類似正態(tài)分布，在同一流場(chǎng)狀態(tài)下，在不同位置釋放的水滴破碎后產(chǎn)生子水滴，其粒徑分布十分近似。

圖7 總收集系數(shù)變化的比例Fig.7 Changing ratio of total collection coefficients

圖8 不同位置釋放水滴的子水滴粒徑分布Fig.8 Diameter distribution of sub-droplet dropping at different positions releasing droplet

圖9給出的是3種不同狀態(tài)下，在遠(yuǎn)場(chǎng)的y=0處釋放的直徑為300 μm的水滴簇破碎后子水滴的粒徑分布。可見，當(dāng)水滴直徑相同時(shí)，在不同翼型和不同流場(chǎng)速度條件下得到的子水滴粒徑分布是十分接近的。所以，水滴發(fā)生破碎之后子水滴的粒徑分布主要與SLD的初始直徑有關(guān)，流場(chǎng)速度、翼型、釋放位置等因素影響十分小。

圖10是SLD破碎后產(chǎn)生子水滴撞擊物面的位置分布，水滴直徑為300 μm，流場(chǎng)速度為150 m/s，釋放點(diǎn)位于NACA0015翼型前某位置。子水滴的位置用碰撞點(diǎn)距離未考慮破碎時(shí)的撞擊點(diǎn)(已考慮變形)的物面距離s來(lái)表示。水滴破碎后，子水滴在空氣動(dòng)力作用下，撞擊位置比不考慮破碎的撞擊位置要遠(yuǎn)離駐點(diǎn)，并散落在一定區(qū)域，圖10 中水滴散落的范圍約為整個(gè)撞擊區(qū)域的0.63%?？梢娝纹扑闀?huì)改變撞擊在機(jī)翼表面的水滴大小及分布，并導(dǎo)致物面局部液態(tài)水的再分配，進(jìn)而影響液態(tài)水在機(jī)翼表面的溢流、傳熱等物理過(guò)程。

圖9 不同流場(chǎng)的子水滴粒徑分布Fig.9 Diameter distribution of sub-droplet in different flow fields

圖10 SLD破碎產(chǎn)生子水滴撞擊位置分布Fig.10 Impingement location distribution of sub-droplet after SLD breakup

3 水滴變形與破碎的綜合影響參數(shù)

通過(guò)前文的分析可知，水滴變形與破碎對(duì)撞擊特性的影響程度隨流速、直徑以及翼型前緣半徑增大而增大。為了綜合考慮各個(gè)影響因素，定義水滴變形與破碎的綜合影響參數(shù)X為

X=v2drL

(19)

式中：rL為機(jī)翼前緣半徑。以X為橫坐標(biāo)，子水滴在翼型表面的覆蓋范圍為縱坐標(biāo)，擬合成直線，如圖11所示。線性回歸的相關(guān)系數(shù)R為0.910 3，說(shuō)明X與子水滴覆蓋范圍有較強(qiáng)的線性相關(guān)關(guān)系。擬合直線的橫截距為0.040 4，未發(fā)生破碎狀態(tài)的X都小于或稍微大于0.040 4，與擬合直線得到結(jié)果吻合。

為了進(jìn)一步驗(yàn)證X作為表征SLD變形與破碎效應(yīng)參數(shù)的適用性，對(duì)X與水滴收集系數(shù)的關(guān)系進(jìn)行了研究。圖12為X與總收集系數(shù)差值的關(guān)系圖。圖中縱坐標(biāo)為考慮變形與破碎和不考慮變形與破碎對(duì)應(yīng)的總收集系數(shù)的差值?？梢钥闯隹偸占禂?shù)的變化量與X呈現(xiàn)較強(qiáng)的線性相關(guān)性，初步說(shuō)明X是適合衡量SLD變形與破碎效應(yīng)的參數(shù)。當(dāng)X=0.040 4時(shí)，考慮變形與破碎和不考慮變形與破碎的總收集系數(shù)差值為0.042，比對(duì)應(yīng)狀態(tài)總收集系數(shù)小一個(gè)數(shù)量級(jí)。水滴直徑為200 μm、流場(chǎng)速度為90 m/s、翼型為NACA0015的狀態(tài)下X=0.040 1，最為接近擬合直線求得的水滴破碎閾值，該狀態(tài)下總收集系數(shù)變化率約為4.56%。說(shuō)明X<0.040 4時(shí)，水滴變形與破碎對(duì)SLD的撞擊特性影響較小。

圖11 X與子水滴覆蓋范圍的關(guān)系Fig.11 Relationship between X and coverage area of sub-droplet

圖12 X與總收集系數(shù)差值的關(guān)系Fig.12 Relationship between X and difference of total collection coefficient

4 結(jié) 論

本文采用數(shù)值計(jì)算的手段研究了水滴直徑、流場(chǎng)速度及繞流物體尺寸變化對(duì)水滴破碎與變形及其效應(yīng)的影響規(guī)律，提出了一個(gè)表征水滴變形與破碎及其效應(yīng)的綜合影響參數(shù)。

1) 水滴直徑、流場(chǎng)速度、機(jī)翼前緣半徑都會(huì)影響水滴的變形與破碎，3個(gè)參數(shù)越大，流場(chǎng)中水滴破碎區(qū)域面積也越大，考慮和未考慮變形與破碎效應(yīng)所對(duì)應(yīng)的水滴收集系數(shù)差異也越大。

2) SLD破碎之后的子水滴分布主要取決于SLD的初始直徑，流場(chǎng)速度、物體尺寸、水滴釋放位置等影響十分小，水滴破碎會(huì)導(dǎo)致物面局部液態(tài)水的再分配，影響液態(tài)水在機(jī)翼表面的溢流、傳熱等物理過(guò)程。

3) 線性回歸分析的結(jié)果說(shuō)明，本文提出的綜合影響參數(shù)X可以體現(xiàn)多因素對(duì)變形與破碎影響的耦合作用，是表征水滴變形與破碎及其效應(yīng)的合適參數(shù)。