基于VOF方法的大水滴袋狀破碎的仿真研究

張文英, 常士楠, 蔣 斌, 雷夢(mèng)龍

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院, 北京 100191)

0 引 言

凍雨的形成是由于大量存在于高空云頂?shù)谋ЯＷ踊蛞簯B(tài)水滴在暖層中凝結(jié)下落，穿過較薄的冷層，由于高層大氣缺少凝結(jié)核而未能發(fā)生相變，仍以液態(tài)水的形式存在。以凍雨或凍毛細(xì)雨形式存在高空中的大水滴(d>50 μm)，我們稱之為過冷大水滴(Supercooled Large Droplets,SLD)。

由于大水滴在空中較之地面的數(shù)密度和雨強(qiáng)貢獻(xiàn)率較大，它對(duì)飛機(jī)表面的撞擊結(jié)冰比常規(guī)尺寸的水滴危害更大[1]。研究常規(guī)尺寸的水滴可以忽略在運(yùn)動(dòng)中的形變和破碎，假設(shè)其始終為球形，且不會(huì)發(fā)生飛濺等現(xiàn)象。而大水滴由于直徑較大，水滴的Oh數(shù)(Ohnesorge數(shù))較小(通常遠(yuǎn)小于0.1)，意味著水滴表面張力和內(nèi)部粘性力的量級(jí)小于所受的氣動(dòng)力，這表征著在水滴撞擊機(jī)翼表面前就會(huì)發(fā)生水滴的變形和破碎，并且在撞擊飛機(jī)迎風(fēng)表面后會(huì)發(fā)生飛濺，生成尺寸更小的子水滴，或?qū)⒃谖丛O(shè)防護(hù)的較大區(qū)域上發(fā)生二次撞擊，出現(xiàn)不可控的結(jié)冰，危害飛機(jī)的安全飛行。由此可見過冷大水滴導(dǎo)致的飛機(jī)結(jié)冰現(xiàn)象較之常規(guī)水滴更為復(fù)雜，引起了各國(guó)研究人員的重視[2-14]。

圍繞運(yùn)動(dòng)液滴的破碎與變形，國(guó)內(nèi)外的研究人員傾注了大量的心血，并取得了一些重要的研究成果。國(guó)外對(duì)于大水滴的實(shí)驗(yàn)研究起步較早，Hinze[2]早在上世紀(jì)中期就開展了在高液-氣密度比和高Re數(shù)條件下液滴破碎的研究，提出了在一定Oh數(shù)條件下，破碎模式隨著We數(shù)改變而變化。Krzeczkowski[3]在Hinze的基礎(chǔ)上進(jìn)行了繼續(xù)研究，確定了破碎模式和We數(shù)的對(duì)應(yīng)關(guān)系，并且分析了典型的袋狀和剪切破碎的形態(tài)特征。Chou和Faeth[4]做了在高液-氣密度比和低Oh數(shù)(Oh?0.1)條件下的液滴破碎的瞬態(tài)特性研究，提出了袋裝破碎時(shí)間與水滴特征時(shí)間相關(guān)的結(jié)論。Rimbert和Castanet[5]認(rèn)為可以用R-T不穩(wěn)定性來解釋袋狀破碎的機(jī)理，或者是由于液滴內(nèi)部的流動(dòng)導(dǎo)致。Guildenbecher[6]對(duì)液滴破碎的實(shí)驗(yàn)文獻(xiàn)做了綜述，并且提出液滴袋狀破碎的機(jī)理是其內(nèi)部液體由兩極向赤道位置的流動(dòng)。

盡管在實(shí)驗(yàn)文章中，研究人員們已經(jīng)給出相對(duì)可靠的破碎模式和We數(shù)對(duì)應(yīng)的范圍，但是由于測(cè)量技術(shù)的限制，難以得到流場(chǎng)中的壓力分布和液滴內(nèi)部的流動(dòng)，因此水滴變形和破碎的機(jī)理尚未明確。伴隨著計(jì)算機(jī)技術(shù)和圖形學(xué)的發(fā)展，研究人員針對(duì)液滴的變形和破碎在數(shù)值仿真方面也取得了一定的成就。Zaleski[7]進(jìn)行了大液滴在Re=1000條件下的數(shù)值仿真，追蹤了相界面的細(xì)節(jié)變化和液滴變形處氣流的渦流結(jié)構(gòu)。Desjardins[8]、Lebas[9]、Sander和Weigand[10]各自做了液滴在不同We數(shù)條件下變形破碎的仿真研究，但是由于研究的重點(diǎn)不同，他們更關(guān)注的是直接數(shù)值仿真或者是大渦模擬的算法應(yīng)用，在防冰以及大水滴的破碎等方面應(yīng)用略顯不足。

國(guó)內(nèi)的相關(guān)研究主要以仿真為主。西北工業(yè)大學(xué)的胡劍平[11]以及上海交通大學(xué)的羅輝[12]等人側(cè)重于應(yīng)用現(xiàn)有的水滴動(dòng)力學(xué)模型去仿真研究SLD對(duì)飛機(jī)部件結(jié)冰的影響。他們的研究結(jié)果證明了，考慮大水滴動(dòng)力學(xué)特性的仿真結(jié)果相較于只考慮常規(guī)尺寸水滴更貼近實(shí)驗(yàn)結(jié)果，充分證明了研究大水滴運(yùn)動(dòng)變形的必要性。然而僅僅從宏觀角度對(duì)水滴運(yùn)動(dòng)行為學(xué)進(jìn)行研究，不能從機(jī)理上解釋水滴變形破碎和飛濺的現(xiàn)象，因此，針對(duì)單個(gè)水滴的運(yùn)動(dòng)和變形破碎進(jìn)行研究是很有意義的。在此方面，北京航空航天大學(xué)的常士楠課題組[13-14]做了關(guān)于水滴袋狀、復(fù)合破碎過程中變形特性研究，對(duì)相關(guān)的理論模型進(jìn)行了改進(jìn)，實(shí)現(xiàn)了水滴變形量和破碎時(shí)間的預(yù)測(cè)。

本文應(yīng)用VOF方法對(duì)不同韋伯?dāng)?shù)下的水滴進(jìn)行數(shù)值模擬，關(guān)注兩相界面變化，得到與實(shí)驗(yàn)一致的變化規(guī)律，應(yīng)用R-T不穩(wěn)定性解釋袋狀破碎現(xiàn)象；統(tǒng)計(jì)了水滴的變形、破碎時(shí)間和運(yùn)動(dòng)速度等瞬態(tài)特性，與前人的實(shí)驗(yàn)數(shù)據(jù)和理論模型加以對(duì)比，揭示其特性與原理。

1 數(shù)值方法與物理模型

在兩相流的數(shù)值仿真方向，主要有兩大類方法：界面追蹤法(Interface tracking Method)和界面捕獲法(Interface capturing Method)，它們分別是基于拉格朗日和歐拉方法。前者將兩相流其中一相視為大量微小粒子的集合，然而它的局限性在于不能自動(dòng)修改拓?fù)浣Y(jié)構(gòu)，并且很難實(shí)現(xiàn)并行運(yùn)算。而基于歐拉方法的界面捕獲法，包括本文采用的VOF(Volume of Fluid)方法和LSM(Level Set Method)方法。VOF方法追蹤每個(gè)網(wǎng)格內(nèi)液體的體積分?jǐn)?shù)，因此有著較好的質(zhì)量守恒特性；但是，在兩相界面陡峭的條件下，LSM方法能夠更好地捕捉界面。但是LSM方法守恒性較差，隨著數(shù)值迭代的積累，在界面變化時(shí)液體體積可能會(huì)丟失，并且它的計(jì)算開銷較之VOF方法更大。綜上所述，本文選擇了VOF方法進(jìn)行仿真。

本文的研究對(duì)象為相對(duì)于靜止空氣有較大的初速度的大水滴，由于與環(huán)境溫差小并且液滴尺寸相對(duì)較小，在研究過程中可以忽略傳熱傳質(zhì)現(xiàn)象，湍流影響和重力作用[15]，因此水滴運(yùn)動(dòng)變形的流場(chǎng)可以視為二維軸對(duì)稱不可壓縮的層流流場(chǎng)，求解相應(yīng)的N-S方程和連續(xù)方程。動(dòng)量方程的控制方程，即考慮表面張力項(xiàng)的N-S方程為：

式中F為液體體積分?jǐn)?shù)，定義如下：

由F定義可知，式(1)中

ρ(F)=ρlF+ρg(1-F)(3)

ρl、ρg分別是液相和氣相的密度。

同理可得

μ(F)=μlF+μg(1-F)(4)

式中μl、μg分別是液相和氣相的動(dòng)力黏度。u為速度矢量，t是時(shí)間，應(yīng)變率張量由D=(u+uT)/2給出。式(1)中最后一項(xiàng)即為表面張力項(xiàng)，σ為表面張力系數(shù)，而κ為氣液界面曲率，δs為表面狄拉克δ函數(shù)。

數(shù)值模擬的計(jì)算模型示意圖如圖1所示，流場(chǎng)為1.5×40 mm的矩形，x軸為軸對(duì)稱條件，其他三個(gè)邊界均為常壓壓力出口條件。初始條件為流場(chǎng)空氣靜止且為常溫常壓(20 ℃，一個(gè)大氣壓)，一個(gè)直徑為0.62 mm的初速度為u0的球形水滴出現(xiàn)在流場(chǎng)中，并從距入口1.5 mm處開始運(yùn)動(dòng)。

圖1 水滴在靜止空氣中運(yùn)動(dòng)破碎示意圖Fig.1 Schematic of a high-speed droplet movement and breakup in still air

影響水滴變形、破碎的參數(shù)主要有水滴初始直徑d0，水滴-氣流相對(duì)速度u，水滴-氣流表面張力σ，水滴、氣流的密度ρl、ρg，動(dòng)力黏度μl、μg等。進(jìn)行無量綱化后得到無量綱參數(shù)表達(dá)式如下：

We=ρgu2d0/σ(5)

Re=ρgud0/μg(6)

采用上述無量綱參數(shù)，可列出算例條件見表1。

由于VOF模型對(duì)于網(wǎng)格質(zhì)量和密度的依賴性較強(qiáng)，本文選用了三套網(wǎng)格進(jìn)行獨(dú)立性驗(yàn)證。在水滴運(yùn)動(dòng)變形的區(qū)域網(wǎng)格尺寸分別是25 μm，10 μm和5 μm，網(wǎng)格數(shù)目分別是25萬，60萬和120萬。驗(yàn)算表明，當(dāng)網(wǎng)格尺寸較大時(shí)，液相界面較粗糙，不能準(zhǔn)確把握變化迅速的水滴變形破碎的特征。但是當(dāng)網(wǎng)格過密時(shí)，增加了大量的計(jì)算時(shí)間且計(jì)算結(jié)果沒有明顯的改善，綜合比較之后，選擇了最小尺寸為10 μm的網(wǎng)格進(jìn)行計(jì)算。

表1 算例初始條件設(shè)置(Oh=0.00473)Table 1 Summary of the simulation initial conditions (Oh=0.00473)

本文對(duì)壓力耦合方程的求解采用了半隱方法，即SIMPLEC算法，在文獻(xiàn)[19]中對(duì)二維層流流動(dòng)的計(jì)算結(jié)果表明，SIMPLEC算法的收斂特性遠(yuǎn)優(yōu)于SIMPLE算法。壓強(qiáng)計(jì)算采用了PRESTO格式，其他物理量采用QUICK格式離散。

2 計(jì)算結(jié)果與分析

2.1 大水滴典型袋狀破碎形態(tài)特征

在不同的We數(shù)下，液滴的破碎呈現(xiàn)出不同的形態(tài)特征，Guildenbecher[6]對(duì)這些破碎模式做了總結(jié)性的分類，它們分別是振動(dòng)破碎、袋狀破碎、復(fù)合模式破碎、剪切破碎和破壞性破碎?？梢杂杀?看出，這些破碎類型之間并不存在We數(shù)的明確界限，不同的研究人員往往給出了不同的結(jié)果。這主要是因?yàn)槟壳皩?duì)破碎模型的定義只是定性而不是定量的，并且由于實(shí)驗(yàn)條件的不同(激波管、風(fēng)洞或者噴嘴)，都對(duì)水滴破碎產(chǎn)生了較大的影響。

對(duì)于表1所列出的速度范圍，即流體相對(duì)速度較低時(shí)，大部分?jǐn)?shù)值分析選擇了層流流動(dòng)進(jìn)行計(jì)算，此時(shí)對(duì)水滴破碎的分析往往基于不穩(wěn)定性理論。對(duì)這一部分的分析往往考慮以下因素：表面張力作用(R-P不穩(wěn)定性)，水滴被環(huán)境空氣加速/減速(R-T不穩(wěn)定性)，以及在液滴邊緣表現(xiàn)出的剪切不穩(wěn)定性(Kelvin-Helmholtz不穩(wěn)定性，簡(jiǎn)稱K-H不穩(wěn)定性)[20]。最后一項(xiàng)常見于剪切破碎模式下的分析，本文暫不考慮。在對(duì)于非紊流袋狀破碎的研究中，R-T不穩(wěn)定性被廣泛應(yīng)用于解釋實(shí)驗(yàn)結(jié)果[21-23]，本文中也將R-T不穩(wěn)定性應(yīng)用于解釋袋狀破碎初期的現(xiàn)象。

表2 不同的液滴破碎模式We數(shù)范圍劃分(Oh數(shù)遠(yuǎn)小于0.1)Table 2 Different breakup regimes as a function of We(Oh is far less than 0.1)

因此，本文針對(duì)不同研究人員確定的袋狀破碎區(qū)間，選擇了We數(shù)為14～22的水滴，作為典型袋狀破碎算例進(jìn)行形態(tài)特征的分析和研究。如圖2所示，可以觀察到水滴在靜止空氣中以中等速度運(yùn)動(dòng)，并逐漸發(fā)生變形破碎的過程。從中可以明顯地觀察到破碎的典型特征和發(fā)生的時(shí)間，在接下來的章節(jié)中將對(duì)這些特征和各類重要的瞬態(tài)性能進(jìn)行分析和研究。 在水滴破碎過程中通常使用特征時(shí)間tc來對(duì)各種瞬態(tài)性能進(jìn)行無量綱化分析，tc的表達(dá)式如下：

圖2 We=16水滴隨時(shí)間變形破碎過程(水滴運(yùn)動(dòng)方向從左至右)各圖對(duì)應(yīng)無量綱時(shí)間為：t/tc=a) 0, b) 0.33, c) 0.55, d) 0.77, e) 1.22, f) 1.77, g) 2.12, h) 2.19.此處tc=452 μs.Fig.2 Numerical simulation for droplet moving from left to right with We=16 at t/tc=a) 0, b) 0.33, c) 0.55, d) 0.77, e) 1.22, f) 1.77, g) 2.12, h) 2.19. (tc=452 μs)

根據(jù)仿真中觀察到的水滴變形和破碎的形態(tài)特征，本文首次提出，將袋狀破碎分為五個(gè)特征階段,具體如下：

1) 梯形變形階段

圖2(a-c)的過程，水滴由球形開始發(fā)生形變。由于水滴運(yùn)動(dòng)方向從左至右，在迎風(fēng)側(cè)的變化較為迅速，由圖3(b)可見，水滴迎風(fēng)側(cè)靠近駐點(diǎn)的位置由于該處的高壓被擠壓變平。而與之相對(duì)應(yīng)的，水滴在背風(fēng)側(cè)的變化較為緩慢，仍基本保持球體不變，這與早期的實(shí)驗(yàn)所觀察到的現(xiàn)象一致[4]。此后由于慣性力和氣動(dòng)阻力的共同作用，水滴整體逐漸被擠壓變平，在無量綱時(shí)間t/tc=0.55時(shí)，可以觀察到明顯的梯形結(jié)構(gòu)。

2) 延伸發(fā)展階段

圖2(c-d)是水滴縱向延伸并且厚度均勻的階段。由圖3(c)可見，在運(yùn)動(dòng)前方駐點(diǎn)處和液體內(nèi)部壓力作用下，水滴在x方向受力逐漸被擠壓。又由于水滴邊緣處空氣的速度增大，靜壓減小，水滴在y方向不斷延展，最終形成了一個(gè)厚度均勻的圓餅形水滴。這個(gè)由球體逐漸變?yōu)楸馄降倪^程在無量綱時(shí)間t/tc=0.77時(shí)，可以顯著觀察到。

3) 水環(huán)形成階段

圖2(d-f)是水環(huán)結(jié)構(gòu)出現(xiàn)的階段。水滴由厚度均勻的餅狀向中心較薄邊緣較厚的袋狀結(jié)構(gòu)發(fā)展，中部的液體向邊緣移動(dòng)，水環(huán)不斷變厚并且作為支撐結(jié)構(gòu)承力，在t/tc=1.77時(shí)可以明顯觀察到水環(huán)的形成。而水滴中心部位反速度方向延展，形成袋狀薄膜。這主要是由于袋狀凹陷處，氣流在駐點(diǎn)處的高壓區(qū)擠壓液體在水滴內(nèi)部由赤道位置向極點(diǎn)位置運(yùn)動(dòng)導(dǎo)致。并且由圖3(e)可以看出，液滴的表面出現(xiàn)了壓力波動(dòng)，然而這些擾動(dòng)還處在表面，由于袋狀薄膜的厚度未達(dá)到臨界值，這些微小的擾動(dòng)沒有穿透液滴。

4) 袋狀破碎階段

圖2(f-g)是袋狀結(jié)構(gòu)不斷向后延展變薄，最終破碎斷裂與水環(huán)分離的過程。這個(gè)過程可以用R-T不穩(wěn)定性來解釋。R-T不穩(wěn)定現(xiàn)象可以簡(jiǎn)單表述為：當(dāng)兩種不同密度流體的分界面做加速運(yùn)動(dòng)且加速度方向垂直于界面時(shí)，此界面是否穩(wěn)定取決于加速度的方向是否從大密度流體指向小密度流體。由于受到氣動(dòng)阻力，并且水滴所受加速度方向由輕質(zhì)流體(空氣)指向重質(zhì)流體(水)，就會(huì)出現(xiàn)R-T不穩(wěn)定現(xiàn)象，表現(xiàn)為袋狀薄膜上出現(xiàn)壓力的擾動(dòng)，見圖3(f-g)。

這個(gè)過程可以用R-T不穩(wěn)定性來解釋，如文獻(xiàn)中所述[21]，假定傳統(tǒng)的二維不穩(wěn)定性理論可以擴(kuò)展至三維條件，根據(jù)R-T理論，可以證得在液滴表面出現(xiàn)特定波長(zhǎng)的擾動(dòng)時(shí)，擾動(dòng)振幅的增長(zhǎng)速度最快，這一特定波長(zhǎng)見式(9)：

其中f是盤狀液滴在靜止空氣中的加速度，其表達(dá)式如式(10)所示：

這種特定波長(zhǎng)的擾動(dòng)迅速在液滴的迎風(fēng)面上形成并發(fā)展，擾動(dòng)的幅度也逐漸增大，形成孔洞最終穿透液膜并迅速擴(kuò)大.在無量綱時(shí)間t/tc=2.12時(shí)，可以觀察到膜狀結(jié)構(gòu)破裂，與水環(huán)分離。

5) 子液滴回縮階段

圖2(g-h)是袋狀薄膜破碎成更小的子液滴的過程。如圖2(h)所示，可以觀察到中心薄膜上和與水環(huán)相接處，均出現(xiàn)了正弦波的形狀。

按照實(shí)驗(yàn)觀測(cè)和R-P不穩(wěn)定性的假設(shè)[20]，受到恒定加速度作用的黏性流體在運(yùn)動(dòng)過程中會(huì)由于不穩(wěn)定性誘導(dǎo)而破碎。即對(duì)于特定形狀的流體，存在一個(gè)特征頻率(fc)的擾動(dòng)，這個(gè)擾動(dòng)的周期(λc)一定而成長(zhǎng)率最大，在運(yùn)動(dòng)過程中占據(jù)統(tǒng)治地位，導(dǎo)致流體發(fā)生破碎。

在這一假設(shè)條件下對(duì)實(shí)驗(yàn)和仿真結(jié)果進(jìn)行分析：Jain[18]在液滴破碎實(shí)驗(yàn)中觀察到子液滴的粒徑呈雙峰分布，分別是水環(huán)破碎和袋狀薄膜破碎產(chǎn)生的，即袋狀薄膜破碎后產(chǎn)生的子液滴粒徑均一。此外，如圖2(g)所示，本文仿真結(jié)果中同樣觀察到，臨界破碎時(shí)薄膜表面出現(xiàn)的正弦波是等振幅的，與前文實(shí)驗(yàn)觀察結(jié)論結(jié)合，可以由理論計(jì)算預(yù)測(cè)袋狀破碎后子水滴的粒徑分布。

在無量綱時(shí)間t/tc=2.19時(shí)，可見到斷裂后的正弦周期狀液面不斷回縮，在表面張力的作用下呈表面積最小化的趨勢(shì)，形成圓形液滴。此外，雖然袋狀薄膜破碎，但外圍的水環(huán)仍保持了原有形狀不變。

圖3 We=16水滴袋裝破碎壓力云圖(水滴運(yùn)動(dòng)方向從左至右)(tc =452 μs)Fig.3 The pressure nephogram of water drop (We=16) during bag breakup(tc =452 μs)

2.2 大水滴典型袋狀破碎瞬態(tài)參數(shù)分析

2.2.1 水滴形變及破碎時(shí)間分析

水滴袋狀破碎過程中的形變和破碎所需的時(shí)間是值得關(guān)注的重要參數(shù)。在前人的工作中，Gordon[16]建立了簡(jiǎn)易的圓柱體水滴模型用于預(yù)測(cè)水滴破碎時(shí)間，模型假設(shè)液體為剛性，氣動(dòng)力將水滴從柱體空間內(nèi)推出即為水滴破碎；此后又有Rourke&Amsden提出基于彈簧振子系統(tǒng)的TAB模型[17]；本課題組創(chuàng)立了將變形過程等效為橢球形水滴的BTB模型[14]，認(rèn)為在變形過程中，當(dāng)半個(gè)水滴的質(zhì)心在直徑方向上發(fā)生一定的相對(duì)位移后破碎。BTB模型引入修正系數(shù)Cm,將相對(duì)變形量D=d/d0表示成We數(shù)的函數(shù)，如式(11)所示。

(D-1+D5-2D-4)=80CmWe(0

本文對(duì)袋狀破碎進(jìn)行仿真研究，并且將仿真結(jié)果與基于BTB模型的計(jì)算結(jié)果和前人的實(shí)驗(yàn)結(jié)果[3-4,18]進(jìn)行比對(duì)，如圖4所示。其中包括來自Krzeczkowski(We=13.5)、Chou &Faeth (We=15)和Mohit Jain(We=20)。圖4中給出了不同We數(shù)下袋狀破碎水滴隨著時(shí)間形變的曲線(實(shí)驗(yàn)為系列點(diǎn))，并以箭頭標(biāo)出了相應(yīng)的水滴破碎點(diǎn)。其中y軸坐標(biāo)為水滴無量綱形變量，即水滴迎風(fēng)直徑d與初始直徑d0的比值(在水滴破碎后，迎風(fēng)直徑視為水環(huán)的外徑)；x軸坐標(biāo)為無量綱時(shí)間，即實(shí)際時(shí)間t與特征時(shí)間tc之比。圖中水滴We數(shù)為16～22，Oh均小于0.005。

圖4 水滴相對(duì)變形量d/d0與無量綱時(shí)間t/tc曲線圖Fig.4 Drop cross-stream diameter as a function of time during bag breakup

首先需要明確的是，水滴破碎變形過程主要受到以下幾種力的作用：氣動(dòng)力，表面張力和黏性力。其中氣動(dòng)力是發(fā)生變形的主因，表面張力和水滴內(nèi)部的黏性力是減緩和阻止水滴變形的作用力。式(5)中給出的常用參數(shù)正是對(duì)這三個(gè)作用力進(jìn)行了無量綱化：We數(shù)代表氣動(dòng)力比表面張力，Oh數(shù)代表了水滴內(nèi)部黏性力比表面張力。在本文關(guān)注的袋狀破碎大水滴領(lǐng)域，水滴Oh數(shù)遠(yuǎn)小于0.01，We>10，即黏性力遠(yuǎn)小于水滴表面張力的量級(jí)，而氣動(dòng)力又遠(yuǎn)大于表面張力，即在水滴運(yùn)動(dòng)的短暫時(shí)間間隔內(nèi)動(dòng)量來不及傳遞，這意味著水滴將先變形而不是改變速度。因此本文分析的重點(diǎn)放在了氣動(dòng)力和水滴形變上，兼顧表面張力的影響。

由圖中三組曲線(仿真組、實(shí)驗(yàn)組和擬合組)可以看出，在每組中曲線中，隨著We數(shù)的增加，水滴的相對(duì)變形量也不斷增大。此外，在仿真組和實(shí)驗(yàn)組的曲線中可以觀察到，水滴破碎后曲線斜率明顯增大，即水滴殘余的水環(huán)迅速擴(kuò)大。聯(lián)系圖3(g,h)的壓力云圖可以分析得出，水滴袋狀薄膜向后凹陷處，存在一個(gè)迎風(fēng)駐點(diǎn)處的高壓區(qū)，而水滴邊緣處由于氣流速度增大而靜壓減小，導(dǎo)致液滴破碎后外側(cè)水環(huán)在y方向受力，有一個(gè)側(cè)向的加速度，快速膨脹。但隨著水滴的破碎，這個(gè)壓差迅速消失，水環(huán)保持一個(gè)勻速的膨脹速度直到破碎。本文之所以關(guān)注此處水環(huán)的變化，是因?yàn)樵谄淦扑楹笮纬傻乃坞m然數(shù)量較少，但卻是子水滴群中體積分?jǐn)?shù)最大的部分，在子水滴群中占統(tǒng)治地位。

對(duì)于水滴袋狀破碎過程中的變形，在Chou和Faeth的實(shí)驗(yàn)文獻(xiàn)中，不僅給出了實(shí)驗(yàn)數(shù)據(jù)，還提出了相關(guān)的經(jīng)驗(yàn)公式，指出在水滴變形初期過程，相對(duì)變形量與無量綱時(shí)間呈線性正相關(guān)，此后為二次函數(shù)關(guān)系，具體公式如式(12)所示：

在對(duì)本文仿真數(shù)據(jù)進(jìn)行擬合后，發(fā)現(xiàn)在We=16的條件下，結(jié)果與經(jīng)驗(yàn)公式吻合很好，如式(13)，且線性段相關(guān)系數(shù)R2=0.995：

但隨著We數(shù)的增長(zhǎng)，線性段曲線斜率增大，逐漸偏離經(jīng)驗(yàn)公式，分析數(shù)據(jù)后認(rèn)為這一經(jīng)驗(yàn)公式更適用于We較低(We=16)的袋裝破碎。并且從圖中數(shù)組曲線可以看出，對(duì)于相對(duì)形變仿真的結(jié)果在全范圍內(nèi)都與試驗(yàn)數(shù)據(jù)吻合較好；而BTB模型主要關(guān)注的是水滴在袋狀薄膜臨界破碎時(shí)期變形量和破碎時(shí)間，因此在破碎后期與實(shí)驗(yàn)數(shù)據(jù)貼合更好。

在水滴破碎時(shí)間的預(yù)測(cè)上，BTB模型效果很好：雖然比起實(shí)驗(yàn)結(jié)果略微偏大，但最大相對(duì)誤差小于15%，并且在We=20的條件下與實(shí)驗(yàn)結(jié)果(圖5)極為貼近。而仿真數(shù)據(jù)則普遍偏小，對(duì)這一結(jié)果進(jìn)行分析，可以參照水滴破碎試驗(yàn)結(jié)果。考慮到實(shí)際實(shí)驗(yàn)中水滴臨界破碎時(shí)，不斷拉伸的袋狀薄膜厚度不斷減小，曲率半徑隨之減小，表面張力不斷增加，與氣動(dòng)力達(dá)到平衡，保持形態(tài)的穩(wěn)定不發(fā)生破碎。直到變形過程中，液膜曲率半徑與水環(huán)半徑一致不再變化，氣動(dòng)力的擾動(dòng)克服了表面張力，薄膜才發(fā)生破碎。但在仿真中，所使用的網(wǎng)格步長(zhǎng)為10 μm，大于了臨界破碎的水膜厚度，在計(jì)算中無法實(shí)現(xiàn)水膜的急劇變薄的過程，導(dǎo)致較之實(shí)驗(yàn)過早地發(fā)生了破碎。對(duì)此應(yīng)做的改進(jìn)是結(jié)合計(jì)算經(jīng)濟(jì)性，在運(yùn)算中使用自適應(yīng)網(wǎng)格，在薄膜處進(jìn)一步細(xì)化網(wǎng)格，才能更好地表現(xiàn)出水膜破碎的完整形態(tài)。

圖5 水滴在We=20條件下的袋狀破碎試驗(yàn)照片[18](氣流方向從右至左, tc=110 μs)Fig.5 Experimental observations[18] for We=20(tc=110 μs)

此外從水滴破碎點(diǎn)中可以看出，隨著We數(shù)的增長(zhǎng)，三組數(shù)據(jù)都顯示出一致的趨勢(shì)：水滴破碎的時(shí)間減小，而臨界破碎的相對(duì)變形量明顯增加。

2.2.2 水滴速度與阻力系數(shù)

在研究袋狀破碎過程中大水滴的速度和受力時(shí)，常常要考慮形變帶來的影響。顯而易見，隨著水滴在不同階段中形狀的改變和迎風(fēng)直徑的增大，受到的阻力比起初始球形有著明顯的增大，在對(duì)阻力和速度進(jìn)行分析時(shí)，常常會(huì)引入阻力系數(shù)CD，將形變與受力以經(jīng)驗(yàn)系數(shù)的形式聯(lián)系起來，CD的表達(dá)式如式(14)所示，其中u是水滴質(zhì)心的瞬時(shí)速度，a是計(jì)算得到的加速度。

圖6給出了不同條件下水滴主體無量綱速度變化量(u0-u)/u0和無量綱時(shí)間t/tc的關(guān)系，以及阻力系數(shù)CD曲線。其中除本文仿真結(jié)果外，還對(duì)照Chou和Faeth[4]的實(shí)驗(yàn)結(jié)果(We=15)以及Jain[18]的三維數(shù)值實(shí)驗(yàn)結(jié)果，驗(yàn)證仿真正確性。

(a)

(b)

從圖6(a)中可以觀察到，液滴質(zhì)心無量綱速度在運(yùn)動(dòng)初期和末期變化平緩，在t/tc=0.6～1.7之間，發(fā)生較為劇烈的變化，即在水滴呈盤形的階段，速度的相對(duì)改變量約為7%。如圖6(b)所示，阻力系數(shù)CD隨著水滴的扁平化，也就是水滴迎風(fēng)直徑的增長(zhǎng)而不斷增大，并且在迎風(fēng)直徑最大時(shí)達(dá)到峰值(CD～1.8)，此后隨著袋狀薄膜開始形成，迎風(fēng)直徑不再顯著變化，而水滴的背風(fēng)面形成了光滑的過渡，水滴前后壓差減小，CD顯著降低。

此外，在起始階段球形液滴的阻力系數(shù)約為0.58，與理論值和實(shí)驗(yàn)數(shù)據(jù)吻合較好。但在峰值處，根據(jù)Chou & Faeth[4]的阻力系數(shù)擬合曲線，袋狀破碎的阻力系數(shù)最大值約為1.6，而仿真的結(jié)果(CD，max～1.8)比起實(shí)驗(yàn)結(jié)果偏大，但與不可壓縮流(Re<2000)中盤形的阻力系數(shù)理論值1.7相近[24]。

3 結(jié) 論

本文應(yīng)用VOF方法，對(duì)大水滴典型袋狀破碎過程中的特性進(jìn)行了探究，得到的主要結(jié)論如下：

1) 袋狀破碎分為五個(gè)特征階段：梯形變形階段、延伸發(fā)展階段、水環(huán)形成階段、袋狀破碎階段和子液滴形成階段；

2) 水膜的破碎可以用R-T不穩(wěn)定性來解釋，子液滴的粒徑可以應(yīng)用R-P不穩(wěn)定性進(jìn)行預(yù)測(cè)；

3) 水滴破碎后殘余的水環(huán)迅速擴(kuò)大。這是由于水滴袋狀薄膜凹陷處，即迎風(fēng)駐點(diǎn)處存在一個(gè)高壓區(qū)，水環(huán)在y向受力導(dǎo)致的；

4) 水滴變形初期過程，相對(duì)變形量與無量綱時(shí)間呈線性正相關(guān)，此后為二次函數(shù)關(guān)系的關(guān)系式，與經(jīng)驗(yàn)公式吻合很好；

5) 水滴破碎時(shí)間隨著We數(shù)的增長(zhǎng)減小，而臨界破碎的相對(duì)變形量明顯增加。

6) 阻力系數(shù)CD隨著水滴迎風(fēng)直徑的增長(zhǎng)而不斷增大，最大值約為1.8，此后顯著降低。