Oldroyd-B 黏彈性液滴碰撞過程的數(shù)值模擬1)

關(guān)新燕 富慶飛,2) 劉 虎 楊立軍

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

? (北京宇航系統(tǒng)控制研究所,北京 100076)

引言

凝膠推進(jìn)劑因具有安全[1-4]、高比沖[5-6]、推力易調(diào)節(jié)[7]、鈍感與環(huán)保[8-9]等優(yōu)點,在導(dǎo)彈武器系統(tǒng)和快速發(fā)射航天運(yùn)載器動力系統(tǒng)等方面有著廣闊的應(yīng)用前景[10-12].但是復(fù)雜的流變特性使其霧化過程存在一定困難.凝膠推進(jìn)劑一般通過添加聚合物膠凝劑配置而成[13-16],制得的推進(jìn)劑具有黏彈性效應(yīng),在霧化過程中會產(chǎn)生黏彈性液絲和液滴等凝膠碎片.

對于牛頓液滴,在各種撞擊條件下進(jìn)行碰撞實驗,得到了合并、反彈、反向分離和拉伸分離的碰撞結(jié)果[17].對于非牛頓液滴的碰撞也已經(jīng)進(jìn)行了一些實驗,Finotello 等[18]通過液滴碰撞實驗研究了剪切變稀非牛頓流體黃原膠的碰撞行為,盡管其中存在著復(fù)雜的碰撞動力學(xué),仍然基于無量綱韋伯?dāng)?shù)(We)和沖擊因子(B)得到了完整的相圖,觀察到了圓盤狀以及碰撞液滴的振蕩行為,這與黏性能量耗散的增加和拉伸作用有關(guān).對于液滴的二次霧化已經(jīng)有相關(guān)的數(shù)值研究[19-23].Khare和Yang[24]使用VOF 捕捉了非牛頓液滴的變形和破裂,分析了冪律型液滴破裂的基本物理原理.發(fā)現(xiàn)當(dāng)液滴向下減速時,會在流動方向上進(jìn)行拉伸,不同的剪切速率導(dǎo)致液滴不同點的黏度不相等,使拉伸變得不對稱,最終液滴形成凹坑直到破裂,破裂時出現(xiàn)串珠結(jié)構(gòu),并且不平衡力產(chǎn)生的扭矩使液滴旋轉(zhuǎn).G u p t a和Sbragaglia[25]數(shù)值研究了牛頓/黏彈性液滴在黏彈性/牛頓流體基質(zhì)中受到剪切時的變形和破裂.采用Lattice-Boltzmann 模型(LBM)模擬兩種具有可變黏度比(液滴與基質(zhì)的黏度比)的不混溶流體,有限差分法(finite difference,FD)用于黏彈性的建模,并使用FENE-P 本構(gòu)方程解釋了聚合物的動力學(xué).結(jié)果發(fā)現(xiàn)在無窮大基質(zhì)中黏彈性的影響很小,而在基質(zhì)局限的情況下會有較明顯的影響.在De(應(yīng)力松弛時間與觀測時間的比值)較小的情況下,與液滴是黏彈性的情況相比,基質(zhì)的黏彈性改變受限液滴穩(wěn)定性的現(xiàn)象更為明顯.Izbassarov和Muradoglu[26]計算研究了黏彈性對固體表面上的液滴撞擊,擴(kuò)散和反彈的影響,采用FENE-CR 模型.發(fā)現(xiàn)更大的Wi(對于特征長度與速度尺度單一的物理問題,Wi與De相一致)的黏彈性更有利于液滴回彈,并且在低韋伯?dāng)?shù)和高雷諾數(shù)下,黏彈性的影響更為明顯.此外,增加平衡接觸角和雷諾數(shù)、減小韋伯?dāng)?shù)都有助于液滴的回彈.通過改變聚合物黏度研究聚合物濃度的影響,結(jié)果是液滴的回彈隨著聚合物黏度的增加而被抑制,這主要?dú)w因于黏性耗散的增強(qiáng).在邊界層中,界面附近的聚合物會產(chǎn)生拉伸應(yīng)力,從而推動接觸線以增加鋪展速率,這是聚合物液滴比牛頓液滴鋪展更廣的主要原因.

文獻(xiàn)[27-28]對冪律液滴相撞進(jìn)行了數(shù)值模擬.Liu 等[27]通過VOF 方法模擬了冪律型凝膠推進(jìn)劑液滴的正向碰撞,研究了反彈、合并和反向分離液滴的現(xiàn)象,結(jié)果表明液滴的反彈取決于韋伯?dāng)?shù)和流體的黏度.在低韋伯?dāng)?shù)下,液滴的最小中心厚度出現(xiàn)晚于其最大變形,這與在高韋伯?dāng)?shù)下的現(xiàn)象相反.而且,復(fù)雜的流場表明最大剪切速率出現(xiàn)在流動被重新定向和加速的點.Sun 等[29]采用LBM 方法在二維坐標(biāo)系中模擬了非牛頓液滴的碰撞,Carreau-Yasuda (CY)模型用于描述非牛頓流變學(xué)對液滴的合并、分離和內(nèi)部混合的影響,We值的范圍大約在20 到200 之間.結(jié)果表明,剪切稀化黏度對碰撞動力學(xué)的影響非常復(fù)雜,并且強(qiáng)烈依賴于流體流變學(xué),隨著非牛頓效應(yīng)的增加,剪切稀化流體促進(jìn)了碰撞液滴的內(nèi)部混合,變形和分離.對于剪切稠化流體,由于隨著剪切增稠而增加的黏性耗散,液滴變形被顯著抑制,促進(jìn)了永久合并.Focke和Bothe[30]采用直接數(shù)值模擬(direct numerical simulation,DNS)和VOF 方法,對Motzigemba 等[31]的實驗數(shù)據(jù)進(jìn)行了研究,結(jié)果表明對于牛頓流體,可以通過選擇有效黏度產(chǎn)生與剪切稀化流體相同的液滴碰撞動力學(xué).他們還提出了一種液膜穩(wěn)定方法[32],以證明高We值下的偏心碰撞也存在有效黏度,非牛頓液滴碰撞會在碰撞復(fù)合體的中心形成薄的液膜,只有在碰撞薄片保持完好無損時,得到的結(jié)果才可靠.Motzigemba等[31]基于流體體積法,對非牛頓液滴的碰撞進(jìn)行了數(shù)值模擬,發(fā)現(xiàn)剪切稀化流體的碰撞液滴變形比牛頓液體的更大.液滴碰撞達(dá)到最大變形所需的無量綱時間與黏度無關(guān),但與韋伯?dāng)?shù)有關(guān),從得到的速度場可以看出,在所有碰撞過程中拉伸流占主導(dǎo).

現(xiàn)有針對黏彈性液滴碰撞的研究大多將其看作純剪切變稀流體,并以冪律模型等對其剪切變稀的流變性質(zhì)進(jìn)行擬合,較少考慮流體具有的彈性的影響.因此本文采用Oldroyd-B 本構(gòu)描述液滴的黏彈性,對兩個等體積液滴的碰撞過程進(jìn)行了直接數(shù)值模擬,研究了松弛時間、黏度比、韋伯?dāng)?shù)和碰撞因子等參數(shù)對黏彈性液滴碰撞過程的影響規(guī)律,以期為凝膠推進(jìn)劑霧化機(jī)理的理解和霧化性能的提升提供參考數(shù)據(jù).

1 數(shù)值計算方法與模型

1.1 數(shù)值算法

數(shù)值模擬采用Basilik 開源代碼[33],將網(wǎng)格自適應(yīng)與流體體積法、連續(xù)表面力模型(continuum surface force,CSF)相結(jié)合,得以有效追蹤界面的運(yùn)動變形和精確計算表面張力.

Basilisk 使用投影法求解流體控制方程(Navier-Stokes 方程)

式中,u=(u,v,w) 為流體的速度,ρ≡ρ(x,t) 為流體的密度,μ≡μ(x,t) 為動力黏度,D為應(yīng)變張量,定義為Dij≡(?iuj+?jui),狄拉克分布函數(shù) δs表示表面張力項集中在界面上,σ 為表面張力系數(shù),κ和n是曲率,并垂直于界面.

對于兩相流,引入第一相的質(zhì)量分?jǐn)?shù)c(x,t),定義密度和黏度為

式中,ρ1,ρ2和μ1,μ2分別是第一相和第二相的密度和黏度.c? 等于c,或者通過對c應(yīng)用一個平滑的空間濾波器構(gòu)造.

密度的平流方程可以用體積分?jǐn)?shù)的等效平流方程來代替

Basilisk 在空間上使用基于四叉樹(二維,三維時為八叉樹)的正交笛卡爾網(wǎng)格對計算域進(jìn)行離散.

對于兩相黏彈性流動的數(shù)值求解,考慮到聚合物的記憶效應(yīng),在N-S 方程中引入聚合物應(yīng)力張量τp

式中,μs是溶液黏度,聚合物張量 τp通常是構(gòu)象張量A的函數(shù)fs(·)

式中,μp為聚合物黏度,λ 為松弛時間.構(gòu)象張量A可以看作是測量聚合物鏈分子變形的內(nèi)部狀態(tài)變量,假設(shè)它對稱且正定,則有

上對流導(dǎo)數(shù)算符 ? 由下式給出

應(yīng)變函數(shù)fS(A)和松弛函數(shù)fR(A)可以通過本構(gòu)方程推導(dǎo)得到,對于Oldroyd-B 本構(gòu)方程有

求解時使用構(gòu)象張量[34-35]的對數(shù) Ψ=lgA代替黏彈性應(yīng)力張量,來解決黏彈性流體數(shù)值求解時常遇到的高維森伯格數(shù)問題[36],即強(qiáng)彈性流動產(chǎn)生高應(yīng)力和精細(xì)特征的區(qū)域時出現(xiàn)數(shù)值不穩(wěn)定性.

1.2 計算模型

圖1 為黏彈性液滴碰撞的仿真模型,初始時刻兩個液滴的直徑為D0,圓心在x方向相距1.05D0,在y方向相距X,以大小相等、方向相反的速度u進(jìn)行碰撞,碰撞因子B=X/D0.計算域為6D0×6D0的正四邊形,初始速度的方向為x方向,與初始速度垂直的方向為y方向.采用Oldroyd-B 模型表示液滴的黏彈性

圖1 黏彈性液滴碰撞的計算模型Fig.1 Computational model of viscoelastic droplet collision

黏度比定義式為β=μs/μ0,其中μ0為零剪切黏度,它是溶液黏度與聚合物黏度之和(μ0=μs+μp).當(dāng)黏彈性液滴進(jìn)行正撞,即B=0 時.如圖2 所示,t=0.01 ms 為液滴初始運(yùn)動時的網(wǎng)格,t=0.05 ms 為液滴即將碰撞時刻的自適應(yīng)網(wǎng)格.對于兩液滴間隙的捕捉,由于本文沒有使用額外的計算模型,因此會使仿真現(xiàn)象產(chǎn)生一定的誤差,未來如果進(jìn)一步加入液滴氣縫的解析模型[37],可能得到更加準(zhǔn)確的結(jié)果.

圖2 不同時刻下液滴撞擊的自適應(yīng)網(wǎng)格Fig.2 Adaptive mesh of droplet collision at different time

圖3 為黏彈性液滴的正撞過程,此時液滴的初始速度為1 m/s,黏度比β=0.1,松弛時間λ=1 ms,圖4 為發(fā)生正撞后液滴尺寸隨時間的變化,這里采用單個液滴的最大寬度Dmax來表示.可以看出液滴從初始時刻的直徑為1 mm 開始,在碰撞擠壓的過程中擴(kuò)散,Dmax會變大,在將近1 ms 時Dmax達(dá)到最大值,此時初始撞擊動能耗盡,除了克服間隙氣體流動耗散掉的能量,動能部分轉(zhuǎn)化為表面能,部分存儲為彈性能.最后速度方向改變,液滴內(nèi)的流體開始向反方向運(yùn)動,直徑又變小,此處稱這一過程為回縮,注意此處的回縮液滴并沒有分離,后面呈現(xiàn)出小幅度振蕩的過程.將液滴的初始速度減小到0.5 m/s,液滴在碰撞后并沒有融合,而是出現(xiàn)反彈的現(xiàn)象.圖5 為黏彈性液滴碰撞反彈的仿真結(jié)果.由于在兩液滴靠近后,其動能較小導(dǎo)致不足以將縫隙中的高壓氣體排除,也沒有等到足夠的時間讓氣體自行排出縫隙,就在互相接觸之前失去了撞擊的慣性,隨后液滴在表面張力的作用下試圖恢復(fù)圓形表面,此時變形了的部分使液滴最終向反方向彈開.

圖3 黏彈性液滴的碰撞融合過程(u=1 m/s)Fig.3 Collision and coalescence process of viscoelastic droplet(u=1 m/s)

圖4 液滴尺寸隨時間的變化曲線Fig.4 Curve of droplet size over time

圖5 黏彈性液滴的碰撞反彈過程(u=0.5 m/s)Fig.5 Collision and bounce process of viscoelastic droplet(u=0.5 m/s)

2 結(jié)果與分析

2.1 松弛時間對黏彈性液滴正撞過程的影響

圖6 為不同松弛時間下液滴最大寬度隨時間的變化.可以看出松弛時間越大,液滴的變形程度越大.在液滴靠近與碰撞的過程中,λ越大時,Dmax的峰值越大,即液滴被擠壓地越扁.這是由于松弛時間增大時,黏彈性應(yīng)力增大,而黏彈性有利于液滴的擴(kuò)散.在液滴回縮的過程中,λ越大時,Dmax的最小值也越小,即在較大的松弛時間下液滴的回縮程度較大,在λ=4 ms 時,Dmax可以減小到初始直徑D0以下.可見不同松弛時間下的液滴回縮階段差異也很大,此時在擠壓階段存儲的彈性能釋放出來,有助于回縮.在液滴振蕩的過程中,λ越大時,振蕩的幅值更大、相應(yīng)的振蕩周期越長.

圖6 不同松弛時間λ 下液滴最大寬度Dmax 隨時間的變化Fig.6 The curve of droplet maximum width Dmax with time at different relaxation time λ

采用如下公式計算液滴的動能KE

式中,ρ為密度,u和v分別為x方向和y方向的速度.在計算模型中通過計算每個網(wǎng)格的動能,然后將所有網(wǎng)格的動能相加來實現(xiàn)系統(tǒng)動能的計算.

液滴表面能SE的計算公式為

這是表面能的二維類比形式.式中σ為表面張力系數(shù),l為氣液界面的總長度.

根據(jù)式(14)和式(15)計算得到不同松弛時間下液滴碰撞過程中動能和表面能隨時間的變化.圖7為不同松弛時間下,系統(tǒng)動能隨時間的變化曲線.從圖中可以看出在系統(tǒng)的初始動能都相同的情況下,松弛時間越大,動能耗散的速率越慢.之后在不同的松弛時間下,動能的變化曲線有相位差,松弛時間增大時,動能變化曲線的拐點向后延遲,振蕩幅值與周期也越大.動能在一開始系統(tǒng)迅速下降到最小值后,又開始上升,松弛時間越大時動能上升的越多,從局部放大圖可以看到在隨后的幾次振蕩中動能的變化也遵循這一規(guī)律.動能峰值隨著振蕩的次數(shù)遞減,由于松弛時間越大時動能耗散的更快,在圖中λ=0.5 ms 時動能只振蕩1 次,λ=1 ms 時動能振蕩2 次.

圖7 不同松弛時間λ 下動能KE 隨時間的變化Fig.7 Kinetic energy KE variation curve with time at different relaxation time λ

圖8 為不同松弛時間下液滴表面能隨時間的變化曲線.在圖中可以看出液滴在擠壓過程中界面總長度增加,表面能迅速上升,松弛時間越大時,表面能的峰值越大,隨后在液滴回縮過程中表面能快速減小,由于此時初始的兩個液滴已經(jīng)合并為一個大液滴,因此其表面能減小到小于初始時刻表面能,在后續(xù)的液滴振蕩過程中隨著液滴的振蕩,相應(yīng)時刻的表面能也有較小的起伏,對于較大的松弛時間,表面能達(dá)到的平衡值相對更小.

圖8 不同松弛時間λ 下表面能SE 隨時間的變化Fig.8 Surface energy SE variation curve with time at different relaxation time λ

Focke和Bothe[30]進(jìn)行了剪切變稀流體的液滴碰撞仿真研究,發(fā)現(xiàn)由于黏性力只對碰撞最初始階段產(chǎn)生影響,因此剪切變稀流體的液滴碰撞過程可以被具有相等有效黏度的牛頓流體復(fù)現(xiàn).對于不同黏度的牛頓液滴碰撞合并過程,Sun 等[38]的仿真結(jié)果表明,當(dāng)Oh增長時,動能的振蕩頻率隨之增大,而振幅相差較小,并且碰撞過程中動能的耗散速度與第一次振蕩的動能變化曲線在不同Oh下沒有顯著區(qū)別,這與本節(jié)得到的改變松弛時間帶來的變化有所不同.

2.2 黏度比對黏彈性液滴正撞過程的影響

圖9 為不同黏度比β時液滴最大寬度隨時間的變化曲線,可以看出β越大時,液滴的碰撞擠壓速率越慢,液滴被擠壓的程度越小,其達(dá)到的最大寬度值越小.隨著β的增大,液滴變形的速率越慢,振蕩幅度更小,更早地達(dá)到平衡,這說明在相同的總黏度下,當(dāng)聚合物黏度的占比降低,即體系中的黏彈性物質(zhì)較少時,回縮與振蕩的趨勢越小.在β=0.2,0.3,0.4 時Dmax的值在回縮后達(dá)到平衡,液滴不再有明顯的振蕩.

圖9 不同黏度比β 下液滴最大寬度Dmax 隨時間的變化Fig.9 The curve of droplet maximum width Dmax with time at different viscosity ratio β

圖10和圖11 分別為不同黏度比β時動能和表面能隨時間的變化曲線.從圖9可以看到對于β越大的液滴,聚合物占比越小,其動能耗散速率越快,后續(xù)動能的振幅更小.而從圖10可以看出液滴的表面能在β的改變下沒有明顯的變化規(guī)律,但可以看出整體上β更大的液滴具有更大的表面能.

圖10 不同黏度比β 下動能KE 隨時間的變化Fig.10 Kinetic energy KE variation curve with time at different viscosity ratio β

圖11 不同黏度比β 下表面能SE 隨時間的變化Fig.11 Surface energy SE variation curve with time at different viscosity ratio β

2.3 韋伯?dāng)?shù)對黏彈性液滴正撞過程的影響

圖12 為不同韋伯?dāng)?shù)下液滴的碰撞過程.在液滴碰撞擠壓階段,韋伯?dāng)?shù)的增長有助于液滴的擠壓,因此達(dá)到的Dmax峰值更大,但由于不同韋伯?dāng)?shù)下液滴的擠壓速率是相等的,所以大韋伯?dāng)?shù)時Dmax達(dá)到峰值的時間延后.由于隨著韋伯?dāng)?shù)的減小,更多的能量被儲存為表面能,因此在收縮階段,更多的表面能釋放,導(dǎo)致液滴收縮得更快,在圖中可以看到更小的韋伯?dāng)?shù)更早完成收縮.在振蕩階段,相比于小韋伯?dāng)?shù)(We=30,50) 的結(jié)果,大韋伯?dāng)?shù)(We=100,200,1000)時融合體的振幅更小.

圖12 不同We 下液滴最大寬度Dmax 隨時間的變化Fig.12 The curve of droplet maximum width Dmax with time at different Weber number

圖13 為不同韋伯?dāng)?shù)下動能隨時間的變化曲線,在動能下降階段,韋伯?dāng)?shù)越小時耗散速率越大,因此液滴的初始碰撞動能更早耗盡,開始回縮.在振蕩階段,小韋伯?dāng)?shù)時的動能更大,因此液滴的振幅越大.圖14 為表面能的變化曲線,可以看出韋伯?dāng)?shù)更大時,表面能越大,其擠壓和回縮過程的相對幅值也更大.

圖13 不同We 下動能KE 隨時間的變化Fig.13 Kinetic energy KE variation curve with time at different Weber number

圖14 不同We 下表面能SE 隨時間的變化Fig.14 Surface energy SE variation curve with time at different Weber number

2.4 碰撞因子對黏彈性液滴碰撞過程的影響

改變碰撞因子B,液滴偏心碰撞時會得到不同的結(jié)果.如圖15 是不同碰撞因子B下兩液滴的碰撞過程,可以看出當(dāng)增大B時,液滴在碰撞后會發(fā)生拉伸與轉(zhuǎn)動,B越大時拉伸得越長.當(dāng)B=0和B=0.2 時,液滴在碰撞擠壓后融合,并在擠壓拉伸過程后較快地恢復(fù)成球形;當(dāng)B=0.4和B=0.6 時,液滴在擠壓后繼續(xù)向著運(yùn)動方向拉伸,可以拉伸到更遠(yuǎn)地距離;當(dāng)B=0.8 時,液滴在擠壓后并沒有發(fā)生融合,而是繼續(xù)向相反的方向運(yùn)動,最終依然分離為兩個液滴.在偏心碰撞過程中,碰撞慣性力的法向分量用來克服液滴之間的氣膜,而切向分量轉(zhuǎn)換成液滴融合體的旋轉(zhuǎn)運(yùn)動,當(dāng)增大B時,更強(qiáng)的切向力使融合體被拉伸.

圖15 不同碰撞因子B 下液滴的碰撞過程Fig.15 The collision process of droplet under different impact factor B

圖16 為不同B下動能隨時間的變化,從中觀察到B越小時動能耗散的速率越快,并且耗散到更低的動能.當(dāng)B=0～ 0.6 時,液滴在其碰撞過程中耗散了絕大部分的動能,但B=0.8 時,液滴在互相接近的過程中只耗散了不到60%的動能,這說明液滴的融合過程會耗散更多的動能.圖17 為對應(yīng)的表面能變化過程,由于B=0.4和B=0.6 時的液滴融合體在計算時間內(nèi)一直處于旋轉(zhuǎn)拉伸狀態(tài),所以其表面能曲線還在上升過程中.B=0.8 時液滴在擠壓過程中表面能達(dá)到最大值,隨后在分離后的球形化過程中降低,并逐漸降到平衡值.

圖16 不同B 下動能KE 隨時間的變化Fig.16 Kinetic energy KE variation curve with time at different impact factor B

圖17 不同B 下表面能SE 隨時間的變化Fig.17 Surface energy SE variation curve with time at different impact factor B

3 結(jié)論

本文對兩個等體積黏彈性液滴的碰撞過程進(jìn)行了直接數(shù)值模擬,考慮了松弛時間λ、黏度比β、韋伯?dāng)?shù)We和碰撞因子B等參數(shù)對黏彈性液滴碰撞過程的影響,主要得到以下結(jié)論.

(1) 黏彈性液滴在正撞融合過程中經(jīng)歷了擠壓擴(kuò)散、回縮和振蕩的過程.

(2) 增大λ有利于液滴擠壓和回縮,可以增大融合體的振蕩振幅和周期,并且延遲變形過程.

(3) 增大β時液滴動能的耗散速度更快,因此其擠壓速度更慢,并且將會阻礙融合體的振蕩.

(4) 增大We時液滴的擠壓程度更大,并且在較小的We下液滴的振幅更大.

(5) 對于側(cè)撞的黏彈性液滴,隨著B的增大,即增大偏心程度時,液滴出現(xiàn)不同的碰撞結(jié)果:B較小時,液滴碰撞后發(fā)生合并;增大B時,合并后的大液滴將發(fā)生旋轉(zhuǎn)與拉伸;繼續(xù)增大B,在靠近后液滴表面互相擠壓,但隨后繼續(xù)運(yùn)動,發(fā)生分離.其中液滴的合并過程將耗散更多的動能.