過冷大水滴變形與破碎的影響因素

李維浩，易賢，李偉斌，王應宇

中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，綿陽 621000

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

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

1 數值計算方法

1.1 空氣流場

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

(1)

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

1.2 水滴軌跡

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

(2)

式中：xd為水滴位移；ρd為水滴密度；Vd為水滴體積；g為重力加速度；Ad為水滴的迎風面積；CD為阻力系數；ua為當地流場速度；ud為水滴速度。引入相對雷諾數：

(3)

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

(4)

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

1.3 水滴變形與破碎

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

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

描述彈簧-質量-阻尼系統的受迫運動的數學方程可以寫為

(5)

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

Y=x/(Cbr)

(6)

(7)

(8)

(9)

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

(10)

其中：

(11)

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

(12)

式中：

(13)

當位移超過初始球體半徑的一半，即Y>1時，認為水滴破碎發生[10]，破碎后產生的子水滴直徑與母水滴直徑的關系表示為

(14)

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

(15)

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

n=D32+2

(16)

子水滴的數量可由質量守恒求得，且認為子水滴是無變形、無振蕩的，并且具有與母水滴相同的速度。

1.4 水滴阻力

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

(17)

水滴變形會引起阻力增加，阻力系數進行以下修正[10]：

CD=CD,S(1+2.632Y)

(18)

1.5 數值計算流程

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

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

圖2 數值計算流程Fig.2 Workflow of numerical calculation

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

2.1 計算條件

SLD在機翼附近變形與破碎和水滴直徑以及水滴和空氣的相對速度有關，而相對速度主要決定于初始來流速度和物體的尺寸外形。對于NACA四位數字翼型，前緣半徑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流場區域劃分

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

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

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

圖4 SLD結冰流場典型區域劃分Fig.4 Typical zoning of SLD icing flow field

2.3 SLD破碎區域面積

表1～表3分別給出了3種翼型流場在不同速度和水滴直徑條件下的破碎區域面積?？梢钥吹?，在所有計算狀態中，只有幾個速度高、水滴直徑大的狀態存在水滴破碎區域。水滴破碎區域面積受到流場速度、水滴直徑以及翼型前緣半徑的共同影響，速度越高、直徑越大、前緣半徑越大，水滴破碎區域面積越大。

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

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

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

2.4 水滴收集特性

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

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

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

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

2.5 子水滴分布

圖8給出的是破碎后產生的子水滴的粒徑分布，母水滴是在同一流場狀態下(v=150 m/s,d=150 μm, NACA0012)，在遠場的不同縱坐標處釋放的，縱坐標y分別為0和0.05?？梢钥闯?，水滴破碎后，會產生大量的子水滴，其分布類似正態分布，在同一流場狀態下，在不同位置釋放的水滴破碎后產生子水滴，其粒徑分布十分近似。

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

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

圖9給出的是3種不同狀態下，在遠場的y=0處釋放的直徑為300 μm的水滴簇破碎后子水滴的粒徑分布?？梢?，當水滴直徑相同時，在不同翼型和不同流場速度條件下得到的子水滴粒徑分布是十分接近的。所以，水滴發生破碎之后子水滴的粒徑分布主要與SLD的初始直徑有關，流場速度、翼型、釋放位置等因素影響十分小。

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

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

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

3 水滴變形與破碎的綜合影響參數

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

X=v2drL

(19)

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

為了進一步驗證X作為表征SLD變形與破碎效應參數的適用性，對X與水滴收集系數的關系進行了研究。圖12為X與總收集系數差值的關系圖。圖中縱坐標為考慮變形與破碎和不考慮變形與破碎對應的總收集系數的差值?？梢钥闯隹偸占禂档淖兓颗cX呈現較強的線性相關性，初步說明X是適合衡量SLD變形與破碎效應的參數。當X=0.040 4時，考慮變形與破碎和不考慮變形與破碎的總收集系數差值為0.042，比對應狀態總收集系數小一個數量級。水滴直徑為200 μm、流場速度為90 m/s、翼型為NACA0015的狀態下X=0.040 1，最為接近擬合直線求得的水滴破碎閾值，該狀態下總收集系數變化率約為4.56%。說明X<0.040 4時，水滴變形與破碎對SLD的撞擊特性影響較小。

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

圖12 X與總收集系數差值的關系Fig.12 Relationship between X and difference of total collection coefficient

4 結 論

本文采用數值計算的手段研究了水滴直徑、流場速度及繞流物體尺寸變化對水滴破碎與變形及其效應的影響規律，提出了一個表征水滴變形與破碎及其效應的綜合影響參數。

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

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

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