受徑向振蕩激勵的黏彈性液滴穩(wěn)定性分析1)

姚慕偉 富慶飛 ?,,2) 楊立軍 ?,

* (中國航發(fā)湖南動力機(jī)械研究所,湖南株洲 412000)

? (北京航空航天大學(xué)宇航學(xué)院,北京 100191)

** (北京航空航天大學(xué)寧波創(chuàng)新研究院先進(jìn)飛行器與空天動力創(chuàng)新研究中心,浙江寧波 315800)

引言

在很多工業(yè)領(lǐng)域都會發(fā)生液滴的不穩(wěn)定過程,如噴墨打印[1-2]、聲懸浮技術(shù)[3-6]、非接觸式測量液體性質(zhì)[7-11]等.Rayleigh[12]和Kelvin[13]最先研究了無黏液滴受小擾動時界面衰減振蕩的自然頻率.在形狀振蕩的過程中,表面張力的作用為使液滴表面恢復(fù)為平衡位置的恢復(fù)力,使液滴的振蕩不斷減小,最終恢復(fù)為球形.隨后Reid[14],Chandrasekhar[15],Miller 和Scriven[16]考慮了自然振蕩中液滴黏性的作用,發(fā)現(xiàn)振蕩衰減的主要原因是在界面附近邊界層內(nèi)的黏性耗散.Bauer 等[17-18]則關(guān)注了流變特性對液滴小振幅振蕩時頻率的影響,發(fā)現(xiàn)流體的黏彈性可能會顯著影響液滴的振蕩行為.Khismatullin 和Nadim[19]對黏彈性液滴形狀振蕩的特征方程進(jìn)行了線性理論分析,發(fā)現(xiàn)對于黏彈性液滴,即使在很小的彈性下也會使高黏度的液滴發(fā)生形狀振蕩.Apfel 等[20-24]對含表面活性劑液滴的形狀振蕩開展了一系列的研究,發(fā)現(xiàn)由于表面活性劑的不均勻分布所導(dǎo)致的Gibbs 彈性是改變自由振蕩頻率最重要的參數(shù).

相比于液滴受小擾動的自由振蕩,Li 等[25]研究了受任意方向加速度的黏性液滴在高速氣流下的R-T 不穩(wěn)定過程,并在忽略液滴黏度的情況下,通過求解特征值問題確定了在中性穩(wěn)定狀態(tài)的一系列臨界Bond 數(shù).通過對遞推關(guān)系的化簡,發(fā)現(xiàn)在大模式數(shù)l條件下增長率與經(jīng)向模數(shù)m無關(guān).吳清等[26]則在Li 等[25]的基礎(chǔ)上考慮了環(huán)境流體的可壓縮性,并將液滴所受的動態(tài)慣性力視為是隨時間變化的.計算結(jié)果顯示雖然氣體的壓縮性具有促進(jìn)液滴破碎的作用,但其影響卻非常微小.Ebo Adou 和Tuckerman[27]首先研究了受時間周期徑向加速度的黏性液滴穩(wěn)定性問題,并進(jìn)行了相關(guān)數(shù)值模擬[28].他們使用了與Kumar 和Tuckerman[29]類似的方法,將Floquet 理論應(yīng)用于球坐標(biāo),與處理水平液面問題所不同的是引入了球諧函數(shù).發(fā)現(xiàn)黏性的存在使中性曲線的不穩(wěn)定舌尖變得光滑,但其分析中卻并未考慮周圍環(huán)境介質(zhì)密度的影響.Li 等[30]則在Ebo Adou 和Tuckerman[27]的基礎(chǔ)上探究了環(huán)境流體密度對穩(wěn)定性的影響,發(fā)現(xiàn)增加環(huán)境流體的密度會使可能激發(fā)的球模態(tài)變窄.Liu 等[31]開展了液滴處于振動隔板上的試驗(yàn),近似計算了液滴上表面駐波的波長,并與Lang’s方程求得的結(jié)果進(jìn)行對比,驗(yàn)證了在低頻下Lang’s方程對球面Faraday 不穩(wěn)定的適用性.

本工作采用Floquet 理論,將Liu 等[31]所研究的牛頓黏性液滴拓展為黏彈性液滴,考慮其受徑向振蕩的加速度情況下的穩(wěn)定性問題,推導(dǎo)得到了表面波增長率隨各參數(shù)變化的系數(shù)矩陣,討論了黏彈性參數(shù)及振蕩條件對液滴不穩(wěn)定性的影響,豐富了液滴Faraday 不穩(wěn)定的理論研究.

1 控制方程

考慮一個黏彈性的不可壓縮球形液滴在靜止氣流中受徑向振蕩的加速度,物理模型如圖1 所示.其加速度只沿著徑向方向,且可以表示成余弦函數(shù)的形式

圖1 液滴受徑向振蕩示意圖Fig.1 Schematic diagram of a droplet subjected to radial oscillation

其中ω表示振蕩頻率,A0表示振蕩幅值.

液滴的線性化控制方程為

式中,ul=(ulr,ulθ,ulφ)表示液滴的擾動速度,pl表示液滴的擾動壓強(qiáng),τl表示液滴的偏應(yīng)力張量.變形的液滴表面為

其中,η表示擾動的液滴表面與平衡位置的偏離量,R為液滴平衡位置的半徑.

類似地,氣相的控制方程為

因?yàn)楸疚难芯康氖欠桥ｎD黏彈性的液滴,采用Jefferys 模型來描述流體的本構(gòu)方程

其中,D=[?ul+(?ul)T]/2 為液體的應(yīng)變率張量,μ0為零剪切黏度,λ1為 應(yīng)力松弛時間,λ2為應(yīng)變馳豫時間.

將擾動量設(shè)為Floquet 解的Fourier 形式

其中,β+iγ 表示Floquet 指數(shù),β為擾動增長率,n表示展開的Fourier 模態(tài)的階數(shù).將擾動量的形式(8)代入本構(gòu)方程(7)中,可以將黏彈性流體的黏度寫成一個和流變模型參數(shù)相關(guān)的有效動力黏度,即

為確定液滴的運(yùn)動,還需要滿足一些邊界條件.首先,速度的徑向分量必須與液滴表面的變形相容,即運(yùn)動邊界條件

其中X1是積分常數(shù),νeff=μeff/ρl為液滴的有效運(yùn)動黏度.在球坐標(biāo)系下的Laplace 方程的通解形式為

考慮到在r→0 時,φl應(yīng)為一有限值,并將φl寫成Floquet 解的Fourier 展開

而ψl是無散度的矢量場,可以用單變量的標(biāo)量函數(shù)ψl(r)表示成以下形式[15]

因?yàn)榧僭O(shè)氣體是無黏的,因此其速度場滿足勢函數(shù)

其中X2是積分常數(shù).考慮到r→∞時 φg→const.則φg的通解滿足

其中Dn是積分常數(shù),可通過運(yùn)動邊界條件確定.將式(39)代入式(13)中,得到Dn的表達(dá)式為

將式(40)代入式(38)和式(39)中,得到pg的表達(dá)式

注意到在式(42)中,除了第一項外,其余項與球諧函數(shù)的參數(shù)l,m和Fourier 模態(tài)的階數(shù)n無關(guān).因此,若要式(42)對于任何l,m及n等于零恒成立,需要滿足

而 式(43)就是所需要的色散方程.

2 線性穩(wěn)定性分析

將色散方程(43)中的系數(shù)化為無量綱的形式

式中λ和q可分別視為是無量綱的激勵頻率的倒數(shù)和無量綱的受迫加速度,c可以用來表征黏度的影響[30].

當(dāng)時,l? 1 有如下Bessel 函數(shù)的漸進(jìn)展開式隨著時間的發(fā)展,如果在液滴表面的擾動幅值一直保持相同的大小,就稱基本流是中性穩(wěn)定的.此時,增長率 β=0,式(45)退化為

根據(jù)方程(53),取多個q的值,得到一系列對應(yīng)的λ,根據(jù)這些(q,λ)的點(diǎn)集做出的曲線即為中性穩(wěn)定曲線.

2.1 中性穩(wěn)定曲線

因?yàn)楸疚牡难芯繉ο鬄榉桥ｎD液滴,因此這里我們主要考慮流變參數(shù)對穩(wěn)定性的影響,即:零剪切黏度μ0、應(yīng)力松弛時間λ1和應(yīng)變馳豫時間λ2.為盡量貼近黏彈性液滴的實(shí)驗(yàn)工況和物性參數(shù),在這里取一組典型的實(shí)驗(yàn)工況.取液滴的半徑為R=5 ×10?4m,所施加的振蕩頻率為ω=5 × 105rad/s,參照Brenn 和Teichtmeister[32]所用的實(shí)驗(yàn)工質(zhì),取質(zhì)量分?jǐn)?shù)為0.8%的普立清2500 溶液,其密度ρl=1000.9 kg/m3,表面張力系數(shù)σ=0.07555 N/m,零剪切黏度μ0=0.7588 Pa·s,應(yīng)力松弛時間λ1=0.163 s.在文獻(xiàn)中并沒有給出λ2的取值,參考Li 等[33]的做法,取λ2=λ1/10,并令為基本工況.因?yàn)槊鞔_了要將μ0=0.7588 Pa·s 代入計算,在這里c的取值就不同于參考文獻(xiàn)[30] 中直接給定的方式.在這里取球波數(shù)l為20,代入式(53)中,再結(jié)合零剪切黏度、密度等參數(shù)即可確定c的值.

根據(jù)基本工況的參數(shù)中性穩(wěn)定曲線的計算結(jié)果如圖2 所示.在圖2 中,實(shí)線表示諧波不穩(wěn)定區(qū)域,虛線表示亞諧波不穩(wěn)定區(qū)域.在基本工況下,第一亞諧波不穩(wěn)定區(qū)域位于 λ ?q平面的第4 象限.根據(jù) λ和q的定義不可能為負(fù)值,所以認(rèn)為在這種條件下第一亞諧波區(qū)域是穩(wěn)定的.不同于無阻尼項的Mathieu 方程,在這里不穩(wěn)定區(qū)域左側(cè)靠近 λ 軸的部分變得光滑,且隨著特征指標(biāo) γ 的增加,距離 λ 軸越來越遠(yuǎn).圖中紅色實(shí)心圓點(diǎn)代表中性穩(wěn)定曲線上q的值最小的點(diǎn),即其橫坐標(biāo)表示整條曲線上的qmin.根據(jù)q的定義,在波數(shù)l確定的情況下,q可以表征振蕩幅值的大小,因此曲線上q最小的點(diǎn)代表使不穩(wěn)定發(fā)生的最小的受迫加速度幅值.在圖2 所示的3 個不穩(wěn)定區(qū)域中,隨著 γ 的增加,qmin逐漸增大,意味著若要在液滴表面激發(fā)出更高階的不穩(wěn)定需要更大的受迫加速度.因此在下面的討論中,僅研究 第一諧波不穩(wěn)定區(qū)域的特征.

圖2 基本工況條件下的中性穩(wěn)定曲線(ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.758 8 Pa·s,λ1=0.163 s,λ2=0.016 3 s)Fig.2 Neutral stability curve under base case condition (ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.758 8 Pa·s,λ1=0.163 s,λ2=0.016 3 s)

圖3 給出了不同零剪切黏度條件下第一諧波區(qū)域的中性穩(wěn)定曲線.隨著零剪切黏度逐漸增大,對應(yīng)的不穩(wěn)定區(qū)域在縮小,不穩(wěn)定舌尖變得越來越光滑,且qmin值在逐漸增大,這意味著隨著液滴零剪切黏度的逐漸增加,整體的穩(wěn)定性增強(qiáng),激勵幅值要響應(yīng)提高才能使液滴發(fā)生諧波模式的振蕩.圖4 顯示了不同應(yīng)力松弛時間條件下液滴的中性穩(wěn)定曲線.當(dāng)應(yīng)力松弛時間逐漸增大時,第一諧波區(qū)域的qmin值逐漸減小,表現(xiàn)為液滴的不穩(wěn)定性增加.

圖3 不同零剪切黏度下的中性穩(wěn)定曲線(ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,=0.7588 Pa·s,λ1=0.163 s,λ2=0.0163 s)Fig.3 Neutral stability curve under different zero-shear viscosity conditions (ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,=0.7588 Pa·s,λ1=0.163 s,λ2=0.0163 s)

圖4 不同應(yīng)力松弛時間下的中性穩(wěn)定曲線(ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.758 8 Pa·s,=0.163 s,λ2=0.0163 s)Fig.4 Neutral stability curve under different stress relaxation time conditions (ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,=0.163 s,λ2=0.0163 s)

圖5 計算了當(dāng)應(yīng)變馳豫時間改變時的中性穩(wěn)定曲線.當(dāng)λ2逐漸增加時,中性穩(wěn)定曲線向右收縮,且qmin值逐漸增大,其變化與零剪切黏度的變化趨勢相同.雖然在前人的研究中[34],與應(yīng)力松弛時間相比,應(yīng)變馳豫時間對流動穩(wěn)定性的影響不大,但在這里卻顯示改變應(yīng)變馳豫時間獲得了更為明顯的變化,推測是因?yàn)樵诒疚闹?基本情況的λ2取得較小,因此成倍地改變λ2影響較為顯著.

圖5 不同應(yīng)變馳豫時間下的中性穩(wěn)定曲線(ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.758 8 Pa·s,λ1=0.163 s,=0.016 3 s)Fig.5 Neutral stability curve under different deformation retardation time conditions (ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.758 8 Pa·s,λ1=0.163 s,=0.016 3 s)

2.2 增長率曲線

在式(46)和式(47)中,q和λ里均含有球波數(shù)l,如果給定一個l的值,在一種確定的試驗(yàn)條件下,q和λ的值是確定的.因此,將球波數(shù)l固定為一個常數(shù),受迫振蕩頻率ω改變,這導(dǎo)致λ隨之變化,在這里,固定A0的值,使得q隨著ω的變化而變化.首先作出在不同加速度振蕩幅值A(chǔ)0條件下的增長率如圖6 所示.對于一個確定的A0,在一特定的頻率段內(nèi),增長率呈階梯狀,隨著激勵頻率的增加,增長率的值逐漸減小.這相當(dāng)于激勵頻率有一個閾值ωcr,當(dāng)ω>ωcr時,液滴才會不穩(wěn)定,且增長率會迅速增加到最大值;而超過這個閾值再進(jìn)一步增加振蕩頻率,液滴雖然還是不穩(wěn)定的,但其表面波的增長卻在逐漸衰減.

圖6 不同加速度振蕩幅值條件下增長率隨振蕩頻率的變化(ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ1=0.163 s,λ2=0.0163 s,l=20)Fig.6 The variation of growth rate versers oscillation frequency under different acceleration conditions (ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ1=0.163 s,λ2=0.0163 s,l=20)

為更好地理解圖6,這里要引入一個“三次方曲線”的概念:在ρg→0 的假設(shè)下,λ與q可通過消去l聯(lián)系起來,得到

對于一種特定的實(shí)驗(yàn)工況,可通過式(54) 在λ?q平面上得到對應(yīng)的曲線,因λ是與q成三次方關(guān)系的,所以這類曲線就叫三次方曲線[30].

在圖7 中,4 條曲線是根據(jù)式(54)做出的不同振蕩幅值的三次方曲線.實(shí)線表示給定增長率的輪廓線,這類似于圖2 中第一諧波區(qū)域中性穩(wěn)定曲線的作圖過程,只不過是β不再視為零而取一些特定的值,這表示特定增長率的輪廓線.隨著振蕩幅值的增加,虛線與不穩(wěn)定區(qū)域相交的部分也越來越寬,這表示不穩(wěn)定也更容易發(fā)生.以增長率為0.25 的綠色輪廓線為例,4 條三次方曲線與該輪廓線的交點(diǎn)分別用黑色實(shí)心圓點(diǎn)和黑色實(shí)心方點(diǎn)標(biāo)出,方點(diǎn)表示三次方曲線進(jìn)入輪廓線時的交點(diǎn),而圓點(diǎn)表示穿出輪廓線時的交點(diǎn).可發(fā)現(xiàn)無論是圓點(diǎn)還是方點(diǎn),隨著加速度振蕩幅值的增大,黑色點(diǎn)的縱坐標(biāo)都是逐漸減小的.聯(lián)系公式(46),縱坐標(biāo)λ可表征無量綱激勵頻率平方的倒數(shù),縱坐標(biāo)的逐漸減小代表激勵頻率逐漸增大,這也就解釋了圖6 中所呈現(xiàn)的振蕩幅值增大,使液滴發(fā)生不穩(wěn)定的振蕩頻率也相應(yīng)增大.

圖7 不同振蕩幅值的三次方曲線與第一諧波不穩(wěn)定區(qū)域輪廓線相交(ρl=1000.9 kg/m3,ρg=0 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ1=0.163 s,λ2=0.0163 s,l=20)Fig.7 The intersection of different oscillation amplitude cubic curves with the contours of the first harmonic instability region (ρl=1000.9 kg/m3,ρg=0 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ1=0.163 s,λ2=0.0163 s,l=20)

圖8 顯示了球波數(shù)l對液滴增長率的影響,在這里及以下的計算中,設(shè)置A0=30 m/s2.隨著球波數(shù)l的增加,液滴的最大增長率減小,且使液滴不穩(wěn)定的振蕩頻率升高.在圖9 中,定義=0.758 8 Pa·s,分別減小和增大零剪切黏度的值,觀察其隨著受迫頻率的變化.零剪切黏度增加時,激勵液滴不穩(wěn)定的頻率值ωcr減小,且最大增長率減小,抑制了擾動的增長,這與前面對于中性穩(wěn)定曲線的分析是一致的.

圖8 不同球波數(shù)l 條件下增長率隨振蕩頻率的變化(ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ1=0.163 s,λ2=0.0163 s,A0=30 m/s2)Fig.8 The variation of growth rate versers oscillation frequency under different spherical wavenumber conditions (ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ1=0.163 s,λ2=0.0163 s,A0=30 m/s2)

圖9 不同零剪切黏度條件下增長率隨振蕩頻率的變化(ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,λ1=0.163 s,λ2=0.0163 s,A0=30 m/s2,l=20)Fig.9 The variation of growth rate versers oscillation frequency under different zero-shear viscosity conditions (ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,λ1=0.163 s,λ2=0.0163 s,A0=30 m/s2,l=20)

圖10 體現(xiàn)了應(yīng)力松弛時間增加時,ωcr逐漸增大.λ1表征流體彈性的大小,應(yīng)力松弛時間越大,液滴更容易發(fā)生失穩(wěn),因此體現(xiàn)在圖10 中為λ1大的情況,最大增長率也較大.圖11 顯示了應(yīng)變馳豫時間增加時,ωcr逐漸減小.在這里λ2的變化趨勢與零剪切黏度改變時對增長率的影響類似,這與在分析中性穩(wěn)定曲線時是相似的.這是因?yàn)樵谑?10) 中,μ0和λ2都是位于分子位置的,改變μ0或λ2時對μeff的影響是同方向的,這也就呈現(xiàn)在增長率或中性穩(wěn)定曲線的結(jié)果中.

圖10 不同應(yīng)力松弛條件下增長率隨振蕩頻率的變化(ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ2=0.0163 s,A0=30 m/s2,l=20)Fig.10 The variation of growth rate versers oscillation frequency under different stress relaxation time conditions (ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ2=0.0163 s,A0=30 m/s2,l=20)

圖11 不同應(yīng)變馳豫時間條件下增長率隨振蕩頻率的變化(ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ1=0.163 s,A0=30 m/s2,l=20)Fig.11 The variation of growth rate versers oscillation frequency under different deformation retardation time conditions (ρl=1000.9 kg/m3,ρg=1.225 kg/m3,σ=0.07555 N/m,μ0=0.7588 Pa·s,λ1=0.163 s,A0=30 m/s2,l=20)

3 結(jié)論

本文通過分析在不同黏彈性參數(shù)下液滴中性穩(wěn)定曲線和增長率的變化趨勢,可以得到以下幾點(diǎn)結(jié)論:

(1) 若要在液滴表面激發(fā)出更高階的不穩(wěn)定,需要更大的受迫加速度;

(2) 零剪切黏度和應(yīng)變馳豫時間的增加,使整體穩(wěn)定性增強(qiáng),若要使液滴發(fā)生諧波模式的不穩(wěn)定需要提高振蕩幅值;應(yīng)力松弛時間的增加則會促進(jìn)液滴的不穩(wěn)定性;

(3) 在一定振蕩幅值情況下,不同振蕩頻率具有不同的擾動增長率,且隨著振蕩頻率的增加,增長率的值逐漸減小;

(4) 與中性穩(wěn)定曲線分析的結(jié)果類似,增加液滴的零剪切黏度和應(yīng)變馳豫時間使液滴的最大增長率減小,同時,增加液滴的球波數(shù)也具有相同的效果;應(yīng)力松弛時間增加,最大增長率增加.