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

中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室，綿陽 621000

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

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

1 數(shù)值計算方法

1.1 空氣流場

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

(1)

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

1.2 水滴軌跡

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

(2)

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

(3)

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

(4)

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

1.3 水滴變形與破碎

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

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

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

(5)

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

Y=x/(Cbr)

(6)

(7)

(8)

(9)

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

(10)

其中：

(11)

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

(12)

式中：

(13)

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

(14)

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

(15)

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

n=D32+2

(16)

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

1.4 水滴阻力

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

(17)

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

CD=CD,S(1+2.632Y)

(18)

1.5 數(shù)值計算流程

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

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

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

2 SLD變形與破碎特性計算分析

2.1 計算條件

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

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

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

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

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

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

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

表1～表3分別給出了3種翼型流場在不同速度和水滴直徑條件下的破碎區(qū)域面積?？梢钥吹?，在所有計算狀態(tài)中，只有幾個速度高、水滴直徑大的狀態(tài)存在水滴破碎區(qū)域。水滴破碎區(qū)域面積受到流場速度、水滴直徑以及翼型前緣半徑的共同影響，速度越高、直徑越大、前緣半徑越大，水滴破碎區(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是有/無變形與破碎模型極限軌跡的對比(速度v=150 m/s，水滴直徑d=300 μm，NACA0015翼型)。圖中可以看出添加變形與破碎模型后撞擊極限變小，這是由于水滴變形引起阻力系數(shù)增大，水滴的跟隨性增強(qiáng)，軌跡的偏折程度會變大。圖6給出的是有/無變形與破碎模型的水滴收集系數(shù)β對比，橫坐標(biāo)是碰撞點離駐點的物面距離s與弦長b的比值。其中圖6(a)顯示的是直徑50 μm的SLD在速度150 m/s的流場中NACA0015翼型上的水滴收集系數(shù)，雖然該流場沒有水滴破碎，但由于SLD在運動中會發(fā)生變形，變形導(dǎo)致水滴阻力系數(shù)增大，水滴運動軌跡與空氣流線一致性更高，從而導(dǎo)致局部水收集系數(shù)和撞擊極限都變小。圖6(b)顯示的是在相同流場下直徑300 μm的SLD的水滴收集系數(shù)，該流場包含有較大的水滴破碎區(qū)域，考慮變形與破碎之后，局部水收集系數(shù)和撞擊極限變化程度比圖6(a)狀態(tài)更大。

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

圖5 極限軌跡對比Fig.5 Comparison of limiting trajectories

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

2.5 子水滴分布

圖8給出的是破碎后產(chǎn)生的子水滴的粒徑分布，母水滴是在同一流場狀態(tài)下(v=150 m/s,d=150 μm, NACA0012)，在遠(yuǎn)場的不同縱坐標(biāo)處釋放的，縱坐標(biāo)y分別為0和0.05?？梢钥闯?，水滴破碎后，會產(chǎn)生大量的子水滴，其分布類似正態(tài)分布，在同一流場狀態(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)場的y=0處釋放的直徑為300 μm的水滴簇破碎后子水滴的粒徑分布。可見，當(dāng)水滴直徑相同時，在不同翼型和不同流場速度條件下得到的子水滴粒徑分布是十分接近的。所以，水滴發(fā)生破碎之后子水滴的粒徑分布主要與SLD的初始直徑有關(guān)，流場速度、翼型、釋放位置等因素影響十分小。

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

圖9 不同流場的子水滴粒徑分布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ù)

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

X=v2drL

(19)

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

為了進(jìn)一步驗證X作為表征SLD變形與破碎效應(yīng)參數(shù)的適用性，對X與水滴收集系數(shù)的關(guān)系進(jìn)行了研究。圖12為X與總收集系數(shù)差值的關(guān)系圖。圖中縱坐標(biāo)為考慮變形與破碎和不考慮變形與破碎對應(yīng)的總收集系數(shù)的差值。可以看出總收集系數(shù)的變化量與X呈現(xiàn)較強(qiáng)的線性相關(guān)性，初步說明X是適合衡量SLD變形與破碎效應(yīng)的參數(shù)。當(dāng)X=0.040 4時，考慮變形與破碎和不考慮變形與破碎的總收集系數(shù)差值為0.042，比對應(yīng)狀態(tài)總收集系數(shù)小一個數(shù)量級。水滴直徑為200 μm、流場速度為90 m/s、翼型為NACA0015的狀態(tài)下X=0.040 1，最為接近擬合直線求得的水滴破碎閾值，該狀態(tài)下總收集系數(shù)變化率約為4.56%。說明X<0.040 4時，水滴變形與破碎對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ù)值計算的手段研究了水滴直徑、流場速度及繞流物體尺寸變化對水滴破碎與變形及其效應(yīng)的影響規(guī)律，提出了一個表征水滴變形與破碎及其效應(yīng)的綜合影響參數(shù)。

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

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

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