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

關新燕 富慶飛,2) 劉 虎 楊立軍

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

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

引言

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

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

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

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

1 數值計算方法與模型

1.1 數值算法

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

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

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

對于兩相流,引入第一相的質量分數c(x,t),定義密度和黏度為

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

密度的平流方程可以用體積分數的等效平流方程來代替

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

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

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

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

上對流導數算符 ? 由下式給出

應變函數fS(A)和松弛函數fR(A)可以通過本構方程推導得到,對于Oldroyd-B 本構方程有

求解時使用構象張量[34-35]的對數 Ψ=lgA代替黏彈性應力張量,來解決黏彈性流體數值求解時常遇到的高維森伯格數問題[36],即強彈性流動產生高應力和精細特征的區域時出現數值不穩定性.

1.2 計算模型

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

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

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

圖2 不同時刻下液滴撞擊的自適應網格Fig.2 Adaptive mesh of droplet collision at different time

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

圖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 結果與分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖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 韋伯數對黏彈性液滴正撞過程的影響

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

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

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

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

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

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

圖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 結論

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

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

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

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

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

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