過冷大水滴粒徑分布的歐拉-拉格朗日混合抽樣算法及對冰型影響

韓涵， 李姚， 印子斐， 孔維梁， 劉洪

(上海交通大學航空航天學院， 上海 200240)

飛機結(jié)冰有極大的危險性，容易引起駕駛員的耦合連鎖反應，導致人-機-環(huán)系統(tǒng)失去穩(wěn)定性，從而引發(fā)飛行風險或?qū)е嘛w行事故。其中， 過冷大水滴(supercooled large droplet，SLD)的環(huán)境給飛行安全帶來極大的威脅，自美國鷹航ATR-72-212飛機事故[1]之后，研究者們逐漸開始關注SLD結(jié)冰現(xiàn)象。近幾年的研究主要集中在SLD的機理研究如壁面性質(zhì)對SLD結(jié)冰的影響[2-5]和SLD結(jié)冰理論建模[6]等，以及SLD動力學特性分析[7- 8]、風洞試驗研究[9]、SLD結(jié)冰不確定性分析[10]等，但目前SLD結(jié)冰理論體系尚不完善，如缺少SLD結(jié)冰云霧環(huán)境對SLD成冰的影響研究。2014年美國聯(lián)邦航空管理局(FAA)發(fā)布了25.140修正案和附錄O，正式提出了SLD環(huán)境的適航符合性驗證要求[11]。與附錄C提出的傳統(tǒng)結(jié)冰環(huán)境不同的是，SLD環(huán)境的粒徑分布為復雜的雙峰分布的形式[12]，該分布特性決定了云霧團中不同粒徑水滴的質(zhì)量分數(shù)不同，特別是液態(tài)水含量(liquid water content，LWC)較少的SLD，其稀疏性決定了它們在空間分布上的強隨機性，同時由于SLD的粒徑較大，故與常規(guī)小粒徑的水滴在動力學效應上有顯著的差異，例如破碎、飛濺和反彈等[12]，且SLD 的結(jié)冰條件會使結(jié)冰速率及水滴收集效率提高，從而使溢流更大，產(chǎn)生過多的流向冰，結(jié)冰范圍最大可超過普通云霧結(jié)冰條件的約 20%的弦長位置[11]，這些特點決定了SLD環(huán)境中的結(jié)冰冰型具有復雜性和隨機性[10,13-14]，這增加了冰型預測和防除冰的難度，而有效地模擬真實的結(jié)冰云霧環(huán)境是準確預測冰型從而優(yōu)化飛機防除冰設計的重要前提，但目前飛機防除冰設計和驗證的基礎手段均無法準確模擬SLD環(huán)境粒徑分布特性。

在結(jié)冰試驗中，結(jié)冰風洞針對SLD 結(jié)冰條件的模擬難度大，模擬條件復雜，故SLD粒徑分布的模擬需求對結(jié)冰風洞中的霧化噴嘴的設計提出了新的要求，目前的主要難點集中在使用兩種具有單集中度分布特性、粒徑范圍不同且能覆蓋SLD結(jié)冰條件中滴譜的粒徑范圍的噴嘴進行混合噴射，通過在穩(wěn)定段和收縮段中的湍流混合，使得云霧在靠近下游的試驗段形成具有“雙峰特性”的顆粒度分布，且要滿足附錄O所規(guī)定的4種滴譜[11]。目前中外成熟的結(jié)冰風洞如美國國家航空航天局(NASA)的Glenn 研究中心結(jié)冰研究風洞(IRT)、意大利航天研究中心(CIRA)的結(jié)冰風洞(IWT)、加拿大國家研究委員會(NRC)的高空模擬結(jié)冰風洞(AIWT)，包括中國空氣動力研究與發(fā)展中心(CARDC)的結(jié)冰風洞也均在起步階段，不能完全覆蓋附錄O規(guī)定的粒徑范圍[11]。故目前冰風洞試驗方法尚不能準確模擬SLD的粒徑分布特性。

數(shù)值模擬的方法難度相對于試驗來說較小，成本較低，可操作空間較大，故一部分研究轉(zhuǎn)向用數(shù)值的方法模擬SLD粒徑分布特性。目前主流的結(jié)冰軟件均不能精確模擬寬粒徑域的分布式粒徑的水滴收集和相應的冰型，通常以中值體積直徑(median volume diameter，MVD)簡化復雜的分布。對于水滴多尺度分布的模擬，目前的研究多停留在對粒徑分布函數(shù)的研究上。Alderliesten[15]研究了Rosin-Rammler分布函數(shù)的參數(shù)與MVD之間的關系；Kollar等[16]研究了不同分布函數(shù)形式對于描述同一水滴粒徑分布的差異。對于粒徑分布對冰型的影響目前研究較少且較淺，朱程香等[17]研究了粒徑分布對水滴撞擊特性和冰型的影響，但其研究的粒徑范圍均在50 μm以下，遠不在附錄O所定義的SLD粒徑范圍內(nèi)。

基于準確模擬SLD環(huán)境粒徑分布特性的必要性，同時又考慮到數(shù)值模擬領域中該研究較匱乏，提出了一種SLD粒徑分布的歐拉-拉格朗日混合抽樣算法。首先基于開源計算流體力學(computational fluid dynamics，CFD)的軟件包OpenFOAM[18]開發(fā)3維結(jié)冰軟件，通過標模風洞測試結(jié)果對單一粒徑條件下的2.5維機翼(2維翼型展向拉伸)的結(jié)冰冰型進行驗證，以測試結(jié)冰模型的準確性。通過在水滴追蹤過程中引入符合真實雙峰分布的抽樣以及水滴位置的隨機性來近似模擬真實飛行過程中SLD環(huán)境中的水滴粒徑分布情況及成冰特點，并在此基礎上研究粒徑分布對冰型特征及冰型空間隨機性的影響，為SLD結(jié)冰的冰型預測和防除冰設計提供一定的理論支撐。

1 數(shù)值計算方法

1.1 翼型繞流計算

由于空氣和冰型、水膜之間存在應力和熱流，翼型外部氣流的速度和溫度是求解結(jié)冰的先決條件；另外，水滴的運動軌跡也受到氣流的影響。將每個結(jié)冰時間步內(nèi)的翼型繞流近似為定常流動，通過調(diào)用OpenFOAM的壓力耦合方程組的半隱式方法(semi implicit method for pressure linked equations，SIMPLE)算法[19]得到定常解。翼型繞流的計算網(wǎng)格為自動生成的笛卡爾-貼體網(wǎng)格[20]，如圖1所示。在每個結(jié)冰時間步均根據(jù)當前冰型重新生成網(wǎng)格并計算空氣流場。湍流模型采用的是剪切應力傳輸(shear stress transfer，SST)模型[21]，近壁面采用的是全y+壁面函數(shù)，湍流普朗特數(shù)為0.9。

圖1 笛卡爾-貼體網(wǎng)格Fig.1 Cartesian-body fitted grids

1.2 基于Rosin-Rammler分布函數(shù)的歐拉-拉格朗日混合抽樣與收集方法

基于Rosin-Rammler分布函數(shù)[22]提出SLD分布式水滴的混合收集方法，從而實現(xiàn)SLD粒徑分布的模擬。根據(jù)SLD滴譜“雙峰分布”特性，50 μm左右兩側(cè)的水滴具有明顯的單峰特性[11]，故以粒徑50 μm作為分界線，對小粒徑水滴群和大粒徑水滴群分別進行追蹤收集。具體方法如下：①用最小二乘法對離散的粒徑和LWC進行Rosin-Rammler分布函數(shù)的擬合；②對于粒徑小于50 μm的水滴，考慮到粒徑差異造成的水滴物理特性上的差異較小，故不通過分布函數(shù)進行重新抽樣，僅用輸入時離散的粒徑和LWC去描述小水滴峰；③對于粒徑大于50 μm的SLD，考慮到每個SLD可能具有破碎、飛濺等特殊的物理特性，個體差異較大，故對定義域大于50 μm的分布函數(shù)段進行大量的粒徑抽樣，從而實現(xiàn)大水滴峰的準確描述；④將兩個“峰”的收集結(jié)果(如質(zhì)量、動量、能量等)進行疊加完成SLD完整粒徑分布的水滴收集。

1.2.1 “雙峰分布”函數(shù)的構造方法

原Rosin-Rammler分布函數(shù)的累計體積頻度分布形式為

(1)

考慮到式(1)只具有單峰特性，故用兩個Rosin-Rammler分布函數(shù)疊加的方式去表現(xiàn)SLD 的雙峰特性，引入顆粒群的質(zhì)量分數(shù)ω，其累積體積頻度可表示為

(2)

對式(2)求導，可得到粒徑的體積頻度函數(shù)為

(3)

1.2.2 SLD粒徑分布的抽樣算法

首先根據(jù)輸入的小水滴的離散分布計算出相應的Vsmall，大水滴的粒徑分布需要滿足Vsmall≤Vbig≤1，采用隨機數(shù)的方式計算出大水滴的每個抽樣粒徑d對應的Vbig，其計算公式為

Vbig=1-(1-Vsmall)rand(0,1)

(4)

式(4)中：rand(0，1)為[0，1]的隨機數(shù)。

再通過式(2)利用二分法可反求出抽樣粒徑的值。

1.2.3 SLD粒徑分布的混合收集方法

主要有兩種方法計算水滴軌跡：歐拉兩相流法[23]和拉格朗日法[24]。歐拉兩相流法將水滴看成空間連續(xù)的介質(zhì)，從而建立了水滴場的概念，通過求解建立在控制體內(nèi)的連續(xù)和動量方程得到水滴速度和容積率。拉格朗日法以單個水滴為研究對象，建立它的運動方程，從而求出其運動軌跡。兩種方法在求解不同問題時有各自的優(yōu)勢。歐拉法在穩(wěn)態(tài)條件下求解分布較集中、統(tǒng)一的顆粒群時比拉格朗日法更有優(yōu)勢。拉格朗日法常用來描述分散的顆粒群的運動狀態(tài)，對于顆粒自身的動力學效應，如SLD的破碎、飛濺的描述更加清晰。故采取小粒徑水滴歐拉法和SLD拉格朗日法相結(jié)合的方法。具體方法如下。

(1)根據(jù)設定的離散粒徑將di<50 μm的部分通過歐拉-歐拉方法求解流場中水滴的連續(xù)和動量方程[25]，求出水滴容積αi和水滴速度Udi，并得到在翼面上的邊界值αi0和Udi0。

(2)通過拉格朗日法對擬合的粒徑分布函數(shù)di>50 μm的部分進行大量的粒徑抽樣，求出噴射水滴的運動軌跡，并得到水膜方程的質(zhì)量源項MLagrangian。

(3)將歐拉和拉格朗日算法下壁面收集的質(zhì)量進行疊加得到最終分布粒徑的質(zhì)量收集Mcollected。

(5)

1.3 水膜流動方法

OpenFOAM采用Bai等[26-27]提出的水膜模型和飛濺模型來模擬未凍結(jié)水在結(jié)冰表面的飛濺、流動與傳熱，該模型基于如下假設：水膜厚度通常不超過1 mm，故水膜在翼面法向上的對流忽略不計，質(zhì)量/動量/能量擴散運輸占主導作用，在翼面切向方向上的對流占主導作用，擴散運輸忽略不計。該模型通過在翼面法向上從翼面到水膜表面對輸運方程的積分獲得求解水膜厚度、流動速度和平均焓的二維方程。

機翼表面的水膜運動控制方程如下。

(6)

SρδU,imp+SρδU,spl

(7)

(8)

式中：U為水膜速度；ρ為水膜密度；δ為水膜厚度；?為微分算子；Sρδ,imp和Sρδ,spl分別為水滴撞擊和飛濺的質(zhì)量源項；pimp和pspl分別為水滴撞擊和飛濺的壓力源項；SρδU,imp和SρδU,spl分別為水滴撞擊和飛濺的動量源項；h為水膜的焓；Sρδh,imp和Sρδh,spl分別為水滴撞擊和飛濺的能量源項。

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

式中：mimp,i、mspl,i分別為水滴撞擊壁面引進的質(zhì)量和水滴飛濺損失的質(zhì)量；vimp,i、vspl,i分別為水滴撞擊壁面和水滴飛濺的瞬時速度；n為壁面的法向；Δt為水膜運輸方程離散后的計算時間步長；himp,i、hspl,i分別為水滴撞擊壁面和水滴飛濺產(chǎn)生的平均焓。

在求解水膜方程時，水膜與空氣流場的數(shù)據(jù)交換過程如圖2所示。由于水膜內(nèi)速度、溫度等變化均為線性變化，故只需要一層網(wǎng)格即可進行離散?？諝饬鲌雠c水膜有共同的邊界條件。首先得到空氣流場中水的質(zhì)量、動量和能量于水膜邊界處的值，傳入水膜得到式(6)～式(8)的源項，再通過求解水膜方程傳出水膜厚度、溫度與速度。

細實線為翼面周圍計算空氣流場、水滴流場的網(wǎng)格； 虛線為計算水膜輸運方程的虛擬網(wǎng)格圖2 水膜與空氣流場數(shù)據(jù)交換示意圖Fig.2 Schematic diagram of data exchange between the film layer and the air flow

1.4 結(jié)冰計算

選取Myers模型[28]預測機翼表面結(jié)冰速率。相比于傳統(tǒng)的Messinger模型[29]，Myers模型的優(yōu)勢在于在計算明冰的結(jié)冰增長時，考慮了水膜和冰層內(nèi)部的溫度梯度變化，可以更準確地捕捉到霜冰到明冰的變化瞬間。

結(jié)冰現(xiàn)象發(fā)生時，質(zhì)量和能量守恒方程可分別表示為

(17)

(18)

式中：t為結(jié)冰的計算時間；B和hw分別為冰高和水膜厚度；ρi和ρw分別為冰和水的密度；β為收集率；V為水滴遠場速度；LWC為液態(tài)水含量；lf為單位質(zhì)量的結(jié)冰潛熱；ki和kw分別為冰和水的導熱系數(shù)；T和θ分別為冰層和水膜內(nèi)部的溫度；z為冰生長高度。

式(17)、式(18)考慮了水膜和冰共同存在的情況，實際求解冰高時，需要根據(jù)能量交換判斷是否會出現(xiàn)水膜，即分為兩種情況討論：霜冰和明冰。

在結(jié)冰的初始階段，過冷水滴撞擊到機翼表面時，由于機翼的表面溫度低于水的凝固點，認為水滴不發(fā)生流動，即結(jié)成冰，此時為霜冰。結(jié)霜冰時，機翼表面的熱流平衡方程為

Qk+Qai+Ql-Qci-Qd-Qs-Qcond=0

(19)

式(19)中：Qci為冰層表面與空氣的對流換熱熱流；Qs為升華熱流；Qk為水滴動能；Ql為結(jié)冰潛熱；Qd為水滴顯熱；Qai為空氣對冰層表面的氣動加熱；Qcond為冰層內(nèi)的導熱熱流。

霜冰結(jié)冰速率方程式為

(20)

式(20)中：ρr為霜冰的密度。

在結(jié)冰的第二階段，當結(jié)冰潛熱不足以把所有的水轉(zhuǎn)化成冰時，機翼表面產(chǎn)生溢流，水和冰同時發(fā)生變化，此時結(jié)冰類型轉(zhuǎn)為明冰。為了使霜冰到明冰的冰高、水膜厚度以及結(jié)冰速率的變化是連續(xù)的，引入關鍵冰高Bg和關鍵時間tg作為判斷霜冰轉(zhuǎn)為明冰的依據(jù)，其表達式分別為

(21)

(22)

式中：Tf和Ts分別為冰的熔點和壁面溫度；Qaf和Qcf分別為空氣對水層的氣動加熱和對流換熱熱流；Qe為蒸發(fā)熱流。

當結(jié)冰類型為明冰時，機翼表面的熱流平衡方程為

Qk+Qaf+Ql+Qin-Qcf-Qd-

Qe-Qcond-Qout=0

(23)

與霜冰相比，式(23)多了對流項Qin和Qout，是水膜的流動產(chǎn)生的能量。水膜的結(jié)冰速率方程為

(24)

式(24)中：ρg為明冰的密度；Ta為空氣溫度。

再將式(24)求解得到的冰高代入式(17)積分可得水膜厚度為

(25)

2 單一粒徑分布水滴的冰型校驗

以2.5維平均冰型校驗結(jié)冰模型的準確性。首先將2維翼型展向拉伸一個弦長的長度并對展向網(wǎng)格加密保證足夠的分辨率，再用2維計算工況進行計算得到2.5維冰型。在模擬SLD結(jié)冰時，拉格朗日法以隨機數(shù)種子初始化粒子的發(fā)射位置以滿足水滴的位置隨機性，從而更接近真實冰風洞試驗或飛行器飛行時的環(huán)境特點。計算得到的冰型在展向上會有一定程度的非均勻度，即冰型隨機性，故在展向提取不同截面的冰型，根據(jù)冰高冰角的特征進行平均處理，得到平均冰型，再與試驗冰型進行對比(用歐拉法求解得到的小水滴流場幾乎是均勻的，故冰型隨機性可忽略不計)。計算結(jié)果的驗證是通過和文獻[30]以及文獻[31]中長為0.533 4 m的NACA0012翼型的冰風洞試驗結(jié)果進行對比。

小粒徑水滴計算工況如表1所示，結(jié)果如圖3所示，計算結(jié)果與試驗結(jié)果吻合良好，表明小水滴的收集、溢流傳熱等計算可靠。

SLD計算工況如表2所示，圖4為平均冰型與試驗冰型的對比結(jié)果。算例3的計算結(jié)果與試驗結(jié)果吻合較好，算例4的計算結(jié)果中結(jié)冰厚度比試驗冰型稍大，基本輪廓相似。

表1 小粒徑水滴計算工況Table 1 Working conditions for small droplets

圖3 冰型校驗Fig.3 Ice shape validation

表2 SLD計算工況Table 2 Working conditions for SLDs

圖4 平均冰型校驗Fig.4 Averaged ice shape validation

SLD的冰型在展向有較明顯的隨機性，隨機性結(jié)果如圖5和圖6所示。該隨機性結(jié)果與冰粗糙度相似，如圖7[32]所示，可見模擬結(jié)果與自然結(jié)冰或冰風洞試驗結(jié)果的粗糙度有一定相似。由于目前的結(jié)冰模擬沒有考慮到試驗設備、結(jié)冰風洞等測量誤差，故和真實冰粗糙度還有一定的差別。

圖5 算例3的冰型隨機性Fig.5 Ice shape randomness for case3

圖6 算例4的冰型隨機性Fig.6 Ice shape randomness for case4

圖7 試驗冰的粗糙度[32]Fig.7 Experimental ice roughness[32]

綜上可知，小水滴霜冰和明冰算例的計算結(jié)果均與試驗結(jié)果對比良好，驗證了收集率模型對于小粒徑水滴的適用性和結(jié)冰模型的準確性；SLD算例收集率偏大，導致結(jié)冰量偏多，可能需要針對SLD進行額外調(diào)整水滴的阻力模型。同時，拉格朗日法考慮了水滴在空間的隨機分布特點，計算得到的SLD冰型在展向具有較明顯的隨機性特征，一定程度上模擬了冰粗糙度的形成。

3 粒徑分布特性的模擬及分布對冰型 的影響

3.1 附錄O雙峰分布的擬合與抽樣結(jié)果

美國聯(lián)邦法規(guī)(CFR)第14篇的第25部航空器適航條例的附錄O針對SLD環(huán)境下的適航標準將過冷大水滴結(jié)冰條件分為凍毛毛雨狀態(tài)(freezing drizzle, FZDZ)以及凍雨狀態(tài)(freezing rain，F(xiàn)ZRA)，并給出了4種SLD結(jié)冰條件滴譜。對粒徑和LWC分別進行抽樣并通過式(2)進行累積體積頻度分布函數(shù)的擬合，得到的結(jié)果如圖8所示，可以看出，擬合結(jié)果與SLD結(jié)冰條件滴譜吻合良好，參數(shù)的擬合結(jié)果如表3所示。

圖8 附錄O條件下的粒徑分布函數(shù)擬合Fig.8 Function fitting of droplet size distributions under appendix O conditions

表3 粒徑分布的擬合參數(shù)Table 3 Fitting parameters for droplet size distributions

再對上述4種粒徑分布的分布曲線進行大量的粒徑抽樣，并將水滴噴射到流場中，通過拉格朗日粒子追蹤和統(tǒng)計的方法獲得抽樣粒徑的分布情況，如圖9所示?？梢钥闯?，抽樣的粒徑分布與理論的體積頻度函數(shù)[式(3)]吻合得較好，結(jié)果符合“雙峰分布”的規(guī)律，從圖9(b)、圖9(d)中可以明顯看出“雙峰”的分布形態(tài)，同時可以看出SLD在空間中整體分布較稀疏，體積分數(shù)相對較小，且粒徑域較寬，這種分布規(guī)律符合人們對其物理上的認知。

圖9 粒徑分布的抽樣與分布曲線的對比Fig.9 Comparison between sampling of the droplet size distribution and the distribution curve

3.2 粒徑分布模擬與單/雙MVD模擬的冰型差異

為量化模擬粒徑分布相對單/雙MVD的精度差異，將各粒徑分布計算得到的冰型與對應的單MVD(整體的MVD)和雙MVD(雙峰各自的MVD)條件下算出的冰型進行對比，并分析粒徑分布對冰型特征的影響程度與冰型隨機性規(guī)律。在3.1節(jié)粒徑分布的條件下進行結(jié)冰計算，計算工況如表4所示。

表4 附錄O條件下的計算工況Table 4 Working conditions for Appendix O

圖10 粒徑分布與單/雙MVD平均冰型的對比Fig.10 Averaged ice shape comparison between methods using droplet size distributions and single/double MVD(s)

從圖10可以看出，不同的粒徑分布對冰型特征產(chǎn)生的影響。結(jié)果顯示，相同結(jié)冰條件下，當MVD<40 μm時，冰型對粒徑分布敏感，分布導致的冰型差異大，特別是冰角特征差別顯著；當MVD>40 μm時，冰型對粒徑分布相對不敏感，且MVD越大，粒徑分布導致的冰型差異越小。

圖11 不同粒徑分布下的冰型隨機性Fig.11 Ice shape randomness under various droplet size distributions

圖11給出了冰型在展向空間的隨機性。結(jié)果顯示，相同結(jié)冰條件下，MVD較小的冰型隨機性主要體現(xiàn)在前緣，而MVD較大的冰型隨機范圍較大，翼型兩側(cè)的流向冰區(qū)域也有一定程度的隨機性，同時MVD較大的冰型隨機程度更強。

綜上，冰型特征和冰型隨機性與粒徑分布的MVD和分布方差有緊密聯(lián)系，故利用控制變量方法進一步探究。

3.3 粒徑分布方差對冰型特征與冰型隨機性的影響

圖12 不同方差的粒徑分布抽樣Fig.12 Droplet size distribution sampling with various variances

圖13 不同粒徑分布方差下的平均冰型Fig.13 Averaged ice shapes with various variances of droplet size distributions

圖13給出了不同分布方差下的平均冰型?？梢钥闯觯S著方差的減小，冰角特征逐漸顯著，且在沒有方差的單一粒徑條件下，冰角特征最顯著。分析可得，當粒徑方差減小時，水滴軌跡的差異減小，前緣收集率相對增大，水膜的溢流相對增強，故冰型的冰角相對增大，下冰角體現(xiàn)得較為顯著。

冰型隨機性結(jié)果如圖14所示，可以看出，隨著粒徑分布方差的減小，冰型隨機范圍有隨著溢流的逐漸增強而擴大的趨勢，冰角處體現(xiàn)得最為顯著，而冰型隨機性的程度差異不顯著。

圖14 不同粒徑分布方差下的冰型隨機性Fig.14 Ice shape randomness with various variances of droplet size distributions

3.4 MVD對冰型特征與冰型隨機性的影響

圖15 不同MVD下的平均冰型Fig.15 Averaged ice shapes with various MVDs

保持噴灑質(zhì)量相同，粒徑分布為單一粒徑，改變 MVD。算例7計算工況與算例6一致。圖15給出了不同MVD下的平均冰型。一方面，隨著MVD的增大，整體收集率呈現(xiàn)增大的趨勢，溢流逐漸增強，流向冰區(qū)域逐漸擴大且增厚。另一方面，當MVD增大到一定程度后，冰角不再“生長”，甚至減小，冰角的特征不顯著，且冰厚不再顯著增大。分析可得，當MVD較小時，隨著MVD增大，溢流增強，冰的增量包括冰角的特征幾乎與MVD呈正相關的關系，但MVD增大到一定值后，水的流量大于等于本地的最大結(jié)冰量，且溢流較強，水迅速流到后方，冰角特征被削弱；冰厚由于逐漸達到最大結(jié)冰量故變化較小，但隨著強溢流的影響，流向冰的范圍逐漸擴大。

冰型隨機性結(jié)果如圖16所示，可以看出，隨著MVD的增加，翼型兩側(cè)流向冰區(qū)域的冰型隨機范圍逐漸擴大，且冰型隨機程度也逐漸增強，下翼面體現(xiàn)得最為顯著?？傮w來說，MVD對冰型隨機性影響較大。

圖16 不同MVD下的冰型隨機性Fig.16 Ice shape randomness with various MVDs

3.5 粒徑分布影響的進一步分析

根據(jù)以上研究，定性得出粒徑分布對于冰型特征與冰型隨機性的影響規(guī)律。再結(jié)合附錄O條件下的冰型預測結(jié)果，得出以下結(jié)論。

對于冰型特征，當MVD較小，即小粒徑水滴的質(zhì)量分數(shù)較大時，粒徑分布方差對冰型特征的影響起主導作用。隨著分布方差的增大，水滴軌跡差異變大，收集的質(zhì)量相對減小，溢流減弱，故冰角特征逐漸減弱，有向中間“縮”的趨勢，冰型差異較大，圖10(a)、圖10(c)與圖13呈現(xiàn)出一樣的趨勢，且算例5的結(jié)冰時間和粒徑分布方差都遠大于算例6，故算例5中的冰角差異體現(xiàn)得更顯著，結(jié)論合理；當MVD較大，即SLD的質(zhì)量分數(shù)較大時，MVD對冰型特征的影響起主導作用。隨著MVD的增大，分布方差帶來的收集率差異減小，冰型差異減小，且MVD越大，粒徑分布方差所帶來的冰型特征差異越小，如圖10(b)、圖10(d)所示。

對于冰型隨機性，從圖14和圖16可以看出，冰型隨機區(qū)域隨著溢流增強而增大，故與MVD呈正相關趨勢，與粒徑分布方差呈負相關趨勢。而冰型的隨機程度主要受MVD影響，反而粒徑分布方差對其影響不顯著。從圖11可以看出，MVD越大，冰型隨機程度越強，隨機區(qū)域也越大，與3.4節(jié)中的結(jié)論一致。

4 結(jié)論

提出了一種歐拉-拉格朗日混合水滴抽樣算法，通過數(shù)值模擬的方法實現(xiàn)了寬粒徑域水滴的溢流結(jié)冰的模擬，并研究了SLD粒徑分布對冰型特征和冰型隨機性的影響。得出如下主要結(jié)論。

(1)通過歐拉-拉格朗日混合水滴抽樣算法，將大小水滴的滴譜分別計算，實現(xiàn)了寬粒徑域分布式水滴的追蹤，更好地為SLD冰型預測和防除冰技術奠定基礎。

(2)用拉格朗日法模擬得到的2.5維冰型在展向空間具有一定的隨機性，該隨機性一定程度上模擬了冰粗糙度的形成，可以為氣動分析提供一定的服務。

(3)冰型隨機性與MVD呈正相關趨勢，而與粒徑分布方差呈負相關趨勢。MVD越大，翼面收集到的質(zhì)量越大，溢流也越強，流向冰區(qū)域逐漸擴大，導致冰型隨機區(qū)域擴大，冰型的隨機程度增強；粒徑分布方差越大，水滴軌跡差異越大，收集到的水滴質(zhì)量越少，導致溢流減弱，冰角“縮小”，使得冰型的隨機區(qū)域減小。

(4)與傳統(tǒng)的MVD冰型計算結(jié)果相比，粒徑分布導致的冰型差異較大。單獨考慮粒徑分布方差和MVD的影響：分布方差越大，方差導致的水滴軌跡差異越大，收集的質(zhì)量越小，進而減弱溢流，故冰角有越小的趨勢；MVD越大，溢流越強，故冰角和冰厚先增大，后由于溢流過大，導致冰角略微減小，且流量大于等于本地最大結(jié)冰量之后，冰厚變化減小。綜合考慮兩者的影響：當MVD較小時，粒徑分布方差對冰型特征的影響占主導作用。方差越大，冰型差異越大；當MVD較大時(遠大于50 μm)，MVD對冰型特征的影響占主導作用。MVD越大，分布方差對冰型特征的影響越小，冰型整體差異越小。