張仕環(huán), 龐明軍, 鄭智穎

(1.常州大學(xué) 機(jī)械與軌道交通學(xué)院,江蘇 常州 213164;2.哈爾濱工業(yè)大學(xué) 能源科學(xué)與工程學(xué)院,哈爾濱 150001)

0 引 言

氣液兩相流廣泛存在于自然界和工業(yè)過程中,如血液中的氣泡運(yùn)動[1],石油輸送過程和鼓泡塔內(nèi)的氣泡上浮運(yùn)動等[2]．氣泡的存在會對液相的運(yùn)動、傳熱傳質(zhì)和表觀黏度等特性產(chǎn)生顯著的影響[3]．因此,研究液相中氣泡的運(yùn)動特性對工農(nóng)業(yè)生產(chǎn)和日常生活具有重要意義．工業(yè)過程涉及的液相多為非Newton流體,而黏彈性流體又是一種典型的非Newton流體,因此研究黏彈性流體中的氣泡運(yùn)動特性意義重大．

為了理解氣泡在非Newton流體中的水動力學(xué)特性,前人已經(jīng)開展了大量的研究[4-19],但關(guān)于氣泡在黏彈性流體中水動力學(xué)特性的研究相對較少．特別是國內(nèi)學(xué)者僅開展了為數(shù)不多的實(shí)驗(yàn)研究[20-22],理論和數(shù)值研究鮮見報道．國外學(xué)者Astarita和Apuzzo[23]開創(chuàng)性的實(shí)驗(yàn)發(fā)現(xiàn),在黏彈性流體中,當(dāng)氣泡體積大到一定程度時,其上浮速度會突然增大(所謂的“速度不連續(xù)性階躍”),氣泡形狀呈現(xiàn)為倒置的淚滴狀．Hassager[24]實(shí)驗(yàn)發(fā)現(xiàn)氣泡在黏彈性流體內(nèi)上浮時,其中心尾跡出現(xiàn)了向下流動(即“負(fù)尾跡”)的現(xiàn)象,與鋼球在黏彈性流體中下沉?xí)r的流線形狀相似[25]．此后,一些學(xué)者實(shí)驗(yàn)再現(xiàn)了黏彈性流體中氣泡上浮的典型特征[26-30]．

在數(shù)值研究方面,Wagner等[31]使用格子Boltzmann方法研究了黏彈性流體中二維氣泡的上浮運(yùn)動,再現(xiàn)了實(shí)驗(yàn)中發(fā)現(xiàn)的帶尖尾部．Pillapakkam等[32]使用水平集方法(level set method),研究了三維氣泡在黏彈性O(shè)ldroyd-B液體中的上浮運(yùn)動,再現(xiàn)了上浮速度不連續(xù)階躍的現(xiàn)象．Fraggedakis等[33]基于Phan-Thien Tanner(PTT)黏彈性模型,使用任意Lagrange-Euler法很好地預(yù)測了在達(dá)到臨界體積時氣泡速度突然增加的現(xiàn)象．然而,隨著研究者們對強(qiáng)彈性流體的數(shù)值計算遇到了數(shù)值問題,即高Weissenberg數(shù)問題(the high Weissenberg number problem,HWNP),為了解決該問題眾多學(xué)者進(jìn)行了深入的研究[34-36]．Yuan等[37]在對構(gòu)象張量進(jìn)行對數(shù)化處理后,使用流體體積(VOF)法研究了高Weissenberg數(shù)下氣泡在黏彈性流體中的上浮運(yùn)動特性．Ohta等[38]基于剪切稀化Carreau模型和FENE-CR黏彈性模型考慮流體的剪切稀化和黏彈性特征,使用水平集和流體體積耦合方法(CLVOF)研究了具有剪切稀化特征黏彈性流體中氣泡的上浮運(yùn)動．

綜上所述,盡管國內(nèi)外學(xué)者針對黏彈性流體中氣泡的上浮運(yùn)動特性開展了一定的研究,并得到了一些有價值的結(jié)論,但相關(guān)現(xiàn)象的機(jī)理解釋仍難以定論．再加上黏彈性流體自身的復(fù)雜特性,研究難度相對較大,國內(nèi)關(guān)于這方面的理論和數(shù)值研究鮮見報道．特別是當(dāng)Weissenberg數(shù)太大時,開展黏彈性流體的數(shù)值計算會遇到“高Weissenberg數(shù)問題”,即構(gòu)象張量分量呈指數(shù)增長會導(dǎo)致計算的不穩(wěn)定性．在實(shí)際工業(yè)應(yīng)用中,也經(jīng)常會遇到彈性較弱(Weissenberg數(shù)較小)的黏性流體,比如低分子量(短鏈)或低濃度的聚合物溶液．為了理解氣泡在黏彈性流體中的上浮運(yùn)動特性,使用數(shù)值模擬方法嘗試研究了低Weissenberg數(shù)下松弛時間、黏性力、表面張力和黏度比對氣泡上浮運(yùn)動特性的影響,詳細(xì)分析了氣泡形狀、上浮速度、尾渦和彈性應(yīng)力的分布特征,為后續(xù)深入研究奠定基礎(chǔ)．

1 物理問題和計算模型

1.1 幾何模型

目前公開發(fā)表的數(shù)值研究氣泡在Oldroyd-B黏彈性流體運(yùn)動的文獻(xiàn)中[32,37],氣泡形狀基本是左右對稱的(除了非穩(wěn)態(tài)工況外),且影響氣泡形狀的因素基本位于氣泡輪廓平面內(nèi)．本文通過研究低Weissenberg數(shù)下,黏性力、松弛時間、表面張力和黏度比對氣泡形狀的影響來反應(yīng)氣泡上浮時的運(yùn)動特性,所以將計算域簡化為二維平面時,氣泡的形狀基本不會受到太大影響．盡管對于二維平面假設(shè),圓形氣泡不是真實(shí)的球形氣泡而是一個圓柱體,導(dǎo)致氣泡尾流和上升阻力存在一定的偏差,但總體上二維平面計算能定性地反映一定的流動特征[39-41]．而且目前的計算是在低Weissenberg數(shù)下進(jìn)行的,氣泡上浮速度相對較小,在一定程度上能將流動失真降低到可接受范圍內(nèi)．另外,目前二維的計算結(jié)果,比如氣泡形狀、速度以及流場信息與文獻(xiàn)基本吻合[37]．所以目前二維計算的假設(shè)是合理的,既能極大地減小計算成本,又能較全面地研究不同參數(shù)對氣泡在黏彈性流體中上浮運(yùn)動的影響．

計算區(qū)域設(shè)置為一矩形區(qū)域,如圖1所示．矩形區(qū)域內(nèi)充滿黏彈性流體,初始半徑R=2.5 mm的圓形氣泡受浮力作用從底部自由上?。疄榱讼诿娴挠绊?以及能使氣泡運(yùn)動取得充分發(fā)展,矩形區(qū)域的高和寬分別取H=90R和W=30R,目前計算區(qū)域的尺寸遠(yuǎn)大于文獻(xiàn)[37]的尺寸．氣泡距底壁的距離為h=5R,以消除底部壁面對氣泡初始運(yùn)動的影響．

圖1 計算模型Fig.1 The computation model

1.2 控制方程

目前的計算是在層流、恒溫且氣液兩相均不可壓縮的前提下進(jìn)行的,流體的控制方程有連續(xù)性方程和動量方程,具體如下:

?·u=0,

(1)

(2)

式中u為速度(m/s),p為壓力(Pa),ρ為流體平均密度(kg/m3),t為時間(s),Ds為黏性應(yīng)力張量,Dp為彈性應(yīng)力張量,Fs為表面張力源項(xiàng),g為重力加速度(m/s2)．Ds和Dp的表達(dá)式如下:

(3)

(4)

式中μ為流體平均黏度(kg/(m·s)),μp為聚合物的黏度(kg/(m·s)),Cij為柔性高分子聚合物或表面活性劑分子的構(gòu)象張量,xi為坐標(biāo)(m),λ為聚合物松弛時間(s),δij為Kronecker符號(i=j時值為1,否則值為0),f(r)為Perterlin函數(shù),即

(5)

(6)

式中L為聚合物分子拉伸長度與平衡長度之比,r為聚合物分子的理想長度．

構(gòu)象張量Cij的輸運(yùn)方程為

(7)

(8)

1.3 界面捕捉方法

使用VOF法捕捉氣液相間的界面,VOF法適用于兩種或多種互不相容的流體,通過引入一個或多個體積分?jǐn)?shù)變量來區(qū)分各相,在每個控制體中,所有相的體積分?jǐn)?shù)之和為1．當(dāng)計算區(qū)域的位置與體積分?jǐn)?shù)已知時,流體的熱物性參數(shù)表示為所有相的體積平均值．流場中相分布可以用相函數(shù)F表示,其被定義如下:

(9)

流體平均密度ρ和黏度μ均是氣液相的體積分?jǐn)?shù)平均值,即

ρ=ρL(F)+ρG(1-F),

(10)

μ=μL(F)+μG(1-F),

(11)

式中ρL為液相的密度(kg/m3),ρG為氣相的密度(kg/m3),μL(=μs+μp)為液相的黏度(kg/(m·s)),μG為氣相的黏度(kg/(m·s)),μs為溶劑的黏度(kg/(m·s)),μp為聚合物的黏度(kg/(m·s))．

VOF法通過求解體積分?jǐn)?shù)(或相函數(shù))輸運(yùn)方程來跟蹤界面,方程如下:

(12)

界面采用Youngs[42]提出的分段線性法(PLIC)進(jìn)行重構(gòu),重構(gòu)后的界面不是連續(xù)的．因此為了保證計算精度,對網(wǎng)格大小有嚴(yán)格的要求．

1.4 表面張力模型

表面張力采用Brackbill 等[43]提出的連續(xù)表面力模型(CSF)計算,氣泡表面的力可以用散度定理表示為體積力,將體積力作為源項(xiàng)加入動量方程(2)中,源項(xiàng)Fs表示為

(13)

(14)

1.5 邊界條件和數(shù)值方法

在圖1所示的矩形計算區(qū)域中,上壁面設(shè)置為壓力出口邊界條件,左、右和下壁面均設(shè)置為無滑移壁面．由于結(jié)構(gòu)化網(wǎng)格具有網(wǎng)格質(zhì)量好、數(shù)據(jù)結(jié)構(gòu)簡單且容易實(shí)現(xiàn)區(qū)域邊界擬合等優(yōu)點(diǎn),所以采用結(jié)構(gòu)化網(wǎng)格對整個計算區(qū)域進(jìn)行網(wǎng)格均勻劃分．梯度項(xiàng)使用基于網(wǎng)格中心的最小二乘法進(jìn)行離散;動量方程采用二階迎風(fēng)格式進(jìn)行離散;壓力項(xiàng)采用PRESTO! (pressure staggering option)方法離散;體積分?jǐn)?shù)采用幾何重構(gòu)法求解．時間上的推進(jìn)采用一階隱式格式,壓力和速度場的耦合選用SIMPLE算法．當(dāng)氣泡形狀和上浮速度基本不隨時間變化時,計算視為收斂．所有參數(shù)設(shè)置完成后進(jìn)行混合初始化,然后在計算域中加入氣泡,開始兩相流的計算．

1.6 計算條件和工況設(shè)置

為了研究氣泡在黏彈性流體中的上浮特性,控制液氣密度比為ρL/ρG=1 000和液氣黏度比為μL/μG=100,分別研究了Galileo數(shù)(Ga)、Weissenberg數(shù)(Wi)、Eotvos數(shù)(Eo)和黏度比(β)對氣泡運(yùn)動的影響．Wi表征聚合物松弛時間與對流時間尺度的比值,Ga表征重力與黏性力的比值,Eo表征重力與表面張力的比值,β是溶劑和溶液的黏度比,分別定義如下:

(15)

表1 計算工況Table 1 Design of the computation case

2 結(jié)果分析和討論

2.1 網(wǎng)格尺寸與時間步長的確定

2.2 結(jié)果可靠性驗(yàn)證

為了確保計算的可靠性,本文將參數(shù)為Ga=1.95,Eo=1,Wi=1和β=0.8時的計算結(jié)果與文獻(xiàn)[37]的結(jié)果進(jìn)行了對比．圖4(左側(cè)為當(dāng)前結(jié)果,右側(cè)為文獻(xiàn)[37]結(jié)果)對比了氣泡的形狀,可以看出氣泡形狀高度吻合．圖5對比了氣泡的上浮速度,時間和速度按照文獻(xiàn)[37]進(jìn)行了無量綱化處理(t*=t/(R/Ug)和v*=v/Ug)．從圖中看出,氣泡穩(wěn)定后的上浮速度分布基本一致,加速階段存在輕微的差異,或許是因以下原因造成的:① 本文使用的是二維模型,文獻(xiàn)使用的是三維模型;② 文獻(xiàn)對構(gòu)象張量進(jìn)行了對數(shù)化處理以處理高Weissenberg數(shù)帶來的數(shù)值問題．總而言之,目前的計算方法是可信的,可以開展目前的研究工作．

圖4 氣泡形狀對比 圖5 氣泡上浮速度對比Fig.4 Comparison of bubble shapes Fig.5 Comparison of bubble center velocities

2.3 氣泡變形分析

氣泡形狀會直接影響周圍液相流場的分布,進(jìn)而對相間的傳質(zhì)傳熱效率產(chǎn)生影響．為此,首先分析了各參數(shù)對黏彈性流體中上浮氣泡形狀的影響．圖6給出了黏度比β=0.2時,氣泡上浮運(yùn)動穩(wěn)定后的形狀．對應(yīng)不同的參數(shù),氣泡呈現(xiàn)為圓形、橢圓形、帽形(破碎和未破碎)、倒置淚滴狀以及其他不規(guī)則形狀．隨著Ga數(shù)的增大,氣泡從圓形向帽形轉(zhuǎn)變,并逐漸趨于扁平狀．當(dāng)Ga=2,在低Wi數(shù)和低Eo數(shù)時,氣泡保持為圓形;當(dāng)Eo數(shù)增大時,由于表面張力的減小,氣泡底部受向上射流的影響發(fā)生了輕微的上凹;當(dāng)Wi≥0.5時,氣泡在低Eo數(shù)下,圓形氣泡底部在彈性力的拉伸作用下變長,隨著Eo數(shù)的變大,氣泡下部形成一個帶尖的尾部(所謂的“尾緣尖”),最終呈現(xiàn)為倒置的淚滴狀．當(dāng)Ga=4時,氣泡形狀總體變得扁平(與Ga=2時相比),而且高Wi數(shù)下的“尾緣尖”現(xiàn)象消失;在較大Eo數(shù)下,出現(xiàn)了帽形氣泡．隨著Ga數(shù)的進(jìn)一步增大,當(dāng)Ga=8時,氣泡變得更加扁平;且在Eo=100時,氣泡底部出現(xiàn)了破碎現(xiàn)象,形成了許多衛(wèi)星氣泡．當(dāng)Ga=16,在Eo=1時,出現(xiàn)了非穩(wěn)態(tài)氣泡,氣泡呈不規(guī)則形狀,上浮路徑為“之”字形;隨著Eo數(shù)的增大,氣泡破碎加?。甒i數(shù)的大小表征了液體黏性和彈性的相對強(qiáng)弱．對于目前的研究,當(dāng)Wi=0.1時,液相的彈性相對較弱,氣泡變形與在純黏液體中的變形相似．隨著Wi數(shù)的增大,液相的彈性增強(qiáng),因此氣泡受液相彈性應(yīng)力的影響變得明顯,當(dāng)Ga=2時,氣泡形狀呈現(xiàn)出尾緣尖特征．隨著Eo數(shù)的增大,即表面張力的減小,氣泡保持原有形狀的能力變?nèi)?氣泡形狀從橢圓形向尾部破碎的帽形轉(zhuǎn)變．

(a) Wi=0.1 (b) Wi=0.5 (c) Wi=1圖6 β=0.2時,氣泡形狀隨Ga數(shù),Wi數(shù)和Eo數(shù)的變化Fig.6 Bubble shapes against Ga,Wi and Eo numbers at β=0.2

為了理解黏度比對氣泡形狀的影響,圖7給出了Eo=10時,黏度比β=0.2,0.5,0.8下,氣泡形狀隨Ga數(shù)和Wi數(shù)的變化．結(jié)合圖6可以看出,隨著黏度比β的增大,原本出現(xiàn)“尾緣尖”的工況,氣泡變形為帽形．這是因?yàn)殡S著黏度比β的增大,聚合物的黏度變小,液相的彈性變?nèi)醵ば栽鰪?qiáng),即液相施加在氣泡上的黏性力變大而彈性力變?。詺馀葑冃闻c小Wi數(shù)的工況相似,類似于在純黏流體中的變形．

(a) Wi=0.1 (b) Wi=0.5 (c) Wi=1圖7 Eo=10時,氣泡形狀隨β,Ga數(shù)和Wi數(shù)的變化Fig.7 Bubble shapes against β,Ga and Wi numbers at Eo=10

為了理解氣泡在黏彈性流體中出現(xiàn)的“尾緣尖”現(xiàn)象,給出聚合分子與氣泡作用示意圖(圖8)．氣泡受浮力作用向上運(yùn)動時,由于黏性的影響會拖動氣泡周圍的液體一起向上運(yùn)動,表征液相彈性的聚合物分子會受到拉伸．當(dāng)聚合物分子被拉伸到一定程度時,會恢復(fù)至平衡態(tài)(類似于彈簧),這時液相中會產(chǎn)生一個與氣泡運(yùn)動相反的力,也就是所謂的彈性應(yīng)力．當(dāng)彈性應(yīng)力足夠大時,會通過黏性作用將力施加到氣泡表面上,該力最終作用在氣泡底部正中間,所以該區(qū)域受拉力的作用變形最明顯,形成“尾緣尖”．

圖8 Wi=1和Ga=2時,聚合物分子對氣泡變形的影響Fig.8 Schematic of effects of polymer molecules on bubble deformation at Wi=1 and Ga=2

另外,對于Eo=1和Ga=16所有的工況,氣泡形狀和運(yùn)動均極不規(guī)則(所謂的非穩(wěn)態(tài)氣泡),氣泡沿“之”字形路徑上?。@是因氣泡兩側(cè)尾渦周期性的脫落,在氣泡左右兩側(cè)產(chǎn)生了周期性的升力引起的．文獻(xiàn)[18,43-44]的實(shí)驗(yàn)研究也發(fā)現(xiàn)了類似的現(xiàn)象．

2.4 氣泡上浮速度分析

氣泡的上浮速度(指氣泡內(nèi)所有節(jié)點(diǎn)的平均速度)是研究氣泡水動力學(xué)的一個重要參數(shù),氣泡的上浮速度在一定程度上影響了液相的含氣率,含氣率的高低會影響液相的表觀黏度和相間的傳遞效率[45-46],所以研究氣泡的上浮速度非常必要．為了研究彈性力和表面張力對氣泡上浮速度的影響,圖9給出了Wi數(shù)和Eo數(shù)對氣泡上浮速度的影響．從圖中可以看出,當(dāng)Eo=1時,隨著Wi數(shù)的增大,氣泡上浮速度隨之增大．對應(yīng)圖6給出的氣泡形狀可知,上浮速度增大的原因可能是在彈性應(yīng)力作用下,氣泡下半部分受到拉伸,氣泡尾部變尖,氣泡整體形狀更接近于流線型,邊界層不易分離,氣泡上升過程受到的形狀阻力較?。?dāng)Eo=10和100時,Wi數(shù)對氣泡上浮速度的影響不明顯,這是因?yàn)闅馀菪螤罹尸F(xiàn)為帽形,橫向尺寸相似,所以氣泡受到的形狀阻力也相似,導(dǎo)致上浮速度變化不大．另外,在相同的Wi數(shù)下,隨著Eo數(shù)的增加,氣泡上浮速度先快速減小,然后緩慢減?。@是由于氣泡的形狀先由圓形變?yōu)槊毙巍⒆罱K近似穩(wěn)定為帽形,導(dǎo)致阻力先快速增大然后趨于穩(wěn)定引起的．

(a) Wi=1

(b) Wi=0.5 (c) Wi=0.1圖9 Ga=4和β=0.2時,Wi數(shù)和Eo數(shù)對氣泡上浮速度的影響Fig.9 Effects of Wi and Eo numbers on the bubble rise velocity at Ga=4 and β=0.2

圖10給出Ga數(shù)對氣泡上浮速度的影響(Eo=10,β=0.2)．當(dāng)Wi=0.1和0.5時,隨著Ga數(shù)的增大,穩(wěn)定后氣泡的上浮速度隨之增大．這是由于隨著Ga數(shù)的增大,在重力保持不變的情況下,氣泡受到液相的黏性阻力減小引起的．在較大的Ga數(shù)(Ga=16)下,由于氣泡前一時刻的破碎,導(dǎo)致氣泡上浮速度在穩(wěn)定前出現(xiàn)了一個峰值．當(dāng)Ga數(shù)較小(如Ga=2)時,氣泡上浮速度隨Wi數(shù)的變化表現(xiàn)出一定的差異．當(dāng)Wi=1時,氣泡經(jīng)歷了一個“加速-減速-再加速”的過程,最后氣泡上浮速度超過相同工況下黏度較小(即Ga數(shù)較大)的工況．第一階段氣泡的加速是因浮力作用引起的;第二階段的減速是因氣泡底部受到射流的作用、發(fā)生上凹變形,變形成帽形,氣泡的形狀阻力增大,所以出現(xiàn)了減速現(xiàn)象;第三階段的加速是因液相中聚合物分子的彈性拉伸,氣泡底部受到聚合物恢復(fù)至平衡態(tài)向下的彈性拉力,氣泡從帽形變?yōu)椤暗怪玫臏I滴”狀,形狀阻力又減小,所以氣泡的上浮速度又增大．因氣泡形狀最終穩(wěn)定在阻力最小的狀態(tài),所以氣泡速度高于低黏度工況下的速度．氣泡的變形過程如圖8所示．

(a) Wi=1

(b) Wi=0.5 (c) Wi=0.1圖10 Eo=10和β=0.2時,Ga數(shù)對氣泡上浮速度的影響Fig.10 Effects of the Ga number on the bubble rise velocity at Eo=10 and β=0.2

另外,從圖9可以看出,當(dāng)其他條件相同,Eo=1時,氣泡上浮過程沒出現(xiàn)再加速現(xiàn)象,而是一次性加速達(dá)到穩(wěn)定的峰值;但Eo≥10時,氣泡速度出現(xiàn)了輕微的波動．這說明形狀阻力的變化是導(dǎo)致氣泡出現(xiàn)再次加速的主要原因．Niethammer等[47]報道過類似的現(xiàn)象,他們認(rèn)為第一階段氣泡的加速是因浮力驅(qū)動引起的,這與本文一致,而第二階段氣泡的減速,則被歸因于氣泡表面彈性應(yīng)力層足夠大以至于對氣泡運(yùn)動產(chǎn)生了影響,對于第三階段氣泡的加速他們認(rèn)為當(dāng)氣泡體積大于臨界體積時,氣泡赤道下方的環(huán)向彈性應(yīng)力引起氣泡再次大幅加速[48]．

圖11展示了黏度比β對氣泡上浮速度的影響(Ga=2,Eo=10)．當(dāng)Wi=1時,隨著黏度比的變大,氣泡上浮速度明顯減?。@是因?yàn)轲ざ缺圈碌脑龃髮?dǎo)致黏性力作用增強(qiáng)而彈性力的作用減小,導(dǎo)致氣泡受到的黏性阻力增大,所以上浮速度變?。ㄟ^對比圖中β=0.2所有的工況可以看出,上文中提到的“加速-減速-再加速”現(xiàn)象隨著Wi數(shù)的增大而變得明顯．這足以說明彈性力對氣泡形狀和上升速度的影響非常大．

(a) Wi=1

(b) Wi=0.5 (c) Wi=0.1圖11 Ga=2和Eo=10時,黏度比β對氣泡上浮速度的影響Fig.11 Effects of viscosity ratio β on the bubble rise velocity at Ga=2 and Eo=10

綜合對比圖9—11可以發(fā)現(xiàn),氣泡在黏彈性液體中上浮時,氣泡上浮速度隨時間的變化表現(xiàn)出現(xiàn)“持續(xù)加速到穩(wěn)定”和“加速-減速-再加速到穩(wěn)定”兩種形式,這不同于氣泡在純黏流體中的現(xiàn)象,說明了彈性的影響非常顯著．氣泡在純黏流體中上浮時,普遍會出現(xiàn)“加速-減速-(某些條件下)輕微再加速”,最后取得穩(wěn)定的終端速度[37],這種速度波動是因浮力和阻力之間的動態(tài)平衡引起的．對于黏彈性流體中“持續(xù)加速到穩(wěn)定”現(xiàn)象是因下述原因引起的,即當(dāng)彈性較弱(Wi=0.1)和表面張力較大(Eo≤10)時,氣泡變形輕微,既不會因射流作用尾部迅速上凹,也不會因?yàn)閺椥粤ψ饔贸尸F(xiàn)“尾緣尖”,此時浮力與氣泡阻力之間動態(tài)平衡產(chǎn)生的速度波動現(xiàn)象因彈性力的存在變得不明顯．當(dāng)彈性較強(qiáng)(Wi≥0.5,β=0.2),黏度較大(Ga=2)和表面張力較弱(Eo≥10)時,氣泡底部出現(xiàn)“尾緣尖”導(dǎo)致形狀阻力大幅減小,從而出現(xiàn)了明顯的再次加速．第二階段減速現(xiàn)象仍存在的原因是,當(dāng)彈性聚合物松弛時間較大時,彈性應(yīng)力的釋放存在一定的延遲,氣泡底部受射流的影響仍會存在上凹的歷程．

2.5 氣泡尾渦分析

為了進(jìn)一步分析氣泡在黏彈性流體中的水動力學(xué)特性,下面給出了氣泡的尾渦分布圖．圖中紅線表示氣泡的邊界輪廓．對于所有的工況,氣泡內(nèi)部都形成了兩對對稱的渦,一對較大,另一對較?。疄榱死斫飧鲄?shù)對氣泡尾渦分布的影響,圖12給出了Ga=4和β=0.2時,Wi數(shù)和Eo數(shù)對氣泡尾渦的影響;圖13給出了Eo=10和β=0.2時,Ga數(shù)和Wi數(shù)對氣泡尾渦的影響．

圖12 Ga=4和β=0.2時,Eo數(shù)和Wi數(shù)對氣泡尾流的影響Fig.12 Effects of Eo and Wi numbers on the bubble wake at Ga=4 and β=0.2

圖13 Eo=10和β=0.2時,Ga數(shù)和Wi數(shù)對氣泡尾流的影響Fig.13 Effects of Ga and Wi numbers on the bubble wake at Eo=10 and β=0.2

從圖12中可以看出,當(dāng)Wi數(shù)相同,在Eo=1時,由于氣泡變形較小,氣泡尾部沒有尾渦生成．這是因氣泡形變小(界面曲率大且變化小),在氣泡界面產(chǎn)生的渦量較小,很快被來流帶走,不會聚集在氣泡尾部．當(dāng)Eo=10和100時,氣泡形變明顯、呈現(xiàn)為帽形,界面曲率變化大且流動阻力(主要是形狀阻力)使氣泡上下側(cè)的壓差變大,導(dǎo)致了尾渦的出現(xiàn),且隨著Eo數(shù)的變大,尾渦被逐漸拉長．當(dāng)氣泡形變成帽形時,尾渦的出現(xiàn)導(dǎo)致氣泡的上浮速度大幅減?。畯膱D13可以看出,當(dāng)其他參數(shù)相同時,隨著Ga數(shù)的增大氣泡尾渦明顯增大．隨著Wi數(shù)的增大,氣泡的尾渦尺寸變小,這是因?yàn)閺椥粤υ谝欢ǔ潭壬夏軌蜃柚箽馀菪巫優(yōu)槊毙?減小氣泡的形狀阻力．這也進(jìn)一步解釋了氣泡上浮速度隨著Wi數(shù)的增大而增大的原因．

2.6 氣泡周圍正應(yīng)力分析

因聚合物正應(yīng)力τyy云圖在一定程度上能體現(xiàn)氣泡周圍聚合物濃度和彈性應(yīng)力的作用效果,而且通過觀察正應(yīng)力云圖可以更直觀地理解彈性應(yīng)力的作用規(guī)律．為此,給出了液相的正應(yīng)力τyy分布云圖(圖14)．

(a) Wi=0.1 (b) Wi=0.5 (c) Wi=1圖14 β=0.2時,τyy分布云圖Fig.14 At β=0.2,contour plots of τyy

從圖14可以看出,大正應(yīng)力區(qū)主要集中在氣泡尾部和氣泡赤道位置附近,這些區(qū)域也是剪切較強(qiáng)的區(qū)域,這也再次證明彈性應(yīng)力對氣泡變形的影響不容忽視．當(dāng)其他參數(shù)相同時,隨著Wi數(shù)的增大,正應(yīng)力的大小和作用范圍總體上表現(xiàn)出增大的趨勢．前文提到彈性應(yīng)力使氣泡尾部受到拉伸,氣泡形狀趨于流線形,從而導(dǎo)致上浮速度增大．從正應(yīng)力云圖來看,大彈性應(yīng)力出現(xiàn)在氣泡赤道附近,分析對應(yīng)位置處的流線不難看出,當(dāng)氣泡上浮時,氣泡頂部和周圍的液體向赤道位置匯聚,導(dǎo)致氣泡表面附近存在較大的速度梯度,所以形成較強(qiáng)的剪切作用．受剪切的影響,氣泡赤道附近聚合物分子被向下拉伸,從而在氣泡上產(chǎn)生一個向上的彈性應(yīng)力,也可使氣泡上浮速度變大,如圖8所示．當(dāng)其他參數(shù)相同,改變氣泡的Eo數(shù)時,氣泡的形狀隨之變化,導(dǎo)致正應(yīng)力的大小和位置也發(fā)生變化．當(dāng)Eo數(shù)較小時,氣泡變形較小,大正應(yīng)力區(qū)域主要出現(xiàn)在氣泡兩側(cè)偏上和尾部區(qū)域;當(dāng)Eo數(shù)較大時,氣泡變成帽形,大正應(yīng)力區(qū)域逐漸移動到帽形的裙邊處;當(dāng)氣泡出現(xiàn)破碎時,大正應(yīng)力區(qū)域分布在脫落氣泡的周圍,也出現(xiàn)了大應(yīng)力區(qū)域脫落的現(xiàn)象．對于大Ga數(shù)的工況,正應(yīng)力作用區(qū)域相對變大,但大應(yīng)力區(qū)域減?。环矫媸蔷酆衔镳ざ葴p小的原因;另一方面是小黏性力的作用導(dǎo)致氣泡上浮速度較大,在聚合物分子松弛時間不變的前提下,在聚合物分子單次拉伸恢復(fù)至平衡態(tài)所需的時間內(nèi),氣泡上浮了更遠(yuǎn)的距離,從而導(dǎo)致正(彈性)應(yīng)力作用區(qū)域變大．

3 結(jié) 論

本文基于Oldroyd-B黏彈性模型,運(yùn)用VOF法模擬了氣泡在黏彈性流體中的上浮運(yùn)動．再現(xiàn)了氣泡在黏彈性流體中呈倒置淚滴形狀的現(xiàn)象,解釋了該現(xiàn)象產(chǎn)生原因,并詳細(xì)分析了Ga數(shù)、Eo數(shù)、Wi數(shù)和黏度比β對氣泡形狀和上浮速度的影響．主要結(jié)論總結(jié)如下:

1) 氣泡在黏彈性流體中上浮時,可呈現(xiàn)為圓形、橢圓形、帽形(破碎和未破碎)、倒置淚滴狀和其他不規(guī)則形狀．隨著Ga數(shù)和Eo數(shù)的增大,黏性力和表面張力減小,氣泡變形程度增加;在大Wi數(shù)和小黏度比β下,彈性較強(qiáng),氣泡底部受到強(qiáng)烈的拉伸作用; 在小Wi數(shù)和大黏度比β下,彈性較弱,氣泡形狀與純黏流體中的相似．當(dāng)Ga=16和Eo=1時,氣泡受到周期性升力的作用表現(xiàn)出非穩(wěn)態(tài)形狀; 當(dāng)Ga=2,Wi≥0.5和β=0.2時,氣泡尾部出現(xiàn)了大小不同的“尾緣尖”,該現(xiàn)象是彈性力、黏性力、表面張力和慣性力綜合作用的結(jié)果,而彈性力是“尾緣尖”現(xiàn)象的必要條件．

2) 氣泡在黏彈性液體中的上浮速度普遍大于在純黏液體中的速度．與純黏流體不同,氣泡在黏彈性流體中上浮時,氣泡上浮速度有持續(xù)加速到終端速度和加速-減速-再加速到終端速度兩種形式,后者的再加速現(xiàn)象是因氣泡“尾緣尖”的出現(xiàn)導(dǎo)致氣泡所受形狀阻力減小引起的．

3) 當(dāng)氣泡變成帽形時,因氣泡上下兩側(cè)壓差較大且界面曲率變化較大,在氣泡尾部出現(xiàn)了尾渦,氣泡尾渦隨著Eo數(shù)的增大而拉長．τyy大應(yīng)力區(qū)主要集中在氣泡尾部和赤道處;隨著Wi數(shù)的增大,聚合物松弛時間變長,且氣泡將聚合物分子拉伸長度變大,所以彈性應(yīng)力作用范圍變大．

致謝本文作者衷心感謝常州大學(xué)“SIETP”基金項(xiàng)目(2022-C-06)對本文的資助．