溫度敏感凝膠推進劑中氣泡運動特性研究1)

宋春雨 李 強 周 昊 李夢哿 封 鋒 吳威濤

(南京理工大學機械工程學院,南京 210094)

引言

凝膠推進劑是指將膠凝劑添加到液體推進劑中形成的一種新型推進劑,一般表現出剪切變稀的非牛頓流變特性,且能夠有效地降低液體推進劑的揮發性,是一種能長期保持穩定的推進劑體系[1].剪切稀化流體是最常見的非牛頓流體,其特征是表觀黏度隨著剪切速率的增加而降低,所以可以使用一種考慮了流體零剪切速率和無窮剪切速率下黏度的本構模型即Carreau-Yasuda 模型對凝膠的流動特性進行模擬[2-3].凝膠的剪切稀化特性使其在剪切速率較低時表現出固體特性,在剪切速率較高時表現出流體特性,因而兼具了固體推進劑和液體推進劑的優點[4].本文所研究的凝膠在低剪切速率下并未表現出彈性效應,所以研究中構造的本構模型并未考慮彈性效應.由于凝膠復雜的流變學特性,難以通過簡單的經驗總結方式進行精確控制,因此需要對其流變特性進行系統研究,以指導其制備、貯存、填裝和使用等過程[5].而在凝膠推進劑的制備和加注過程中,氣泡的混入是難以避免的,這可能會導致推力振蕩、推力斷續或熄火等問題.此外,氣泡的混入還會影響推進劑的傳熱性能和泵的氣蝕現象[6].所以對于凝膠推進劑中氣泡運動行為的探究有重要的意義.

大量的學者研究了氣泡在非牛頓流體中的運動情況并探討了流體的流變特性對氣泡運動的影響.Premlata 等[7]的研究發現氣泡在牛頓流體中可以保持直線上升,在非牛頓流體中氣泡的運動軌跡與流體的流變指數與特征時間都有關系.Davidson 等[8]的研究表明,在黏性流體中,氣泡上升的速度與液體的黏度成反比,與液體的密度成正比.同時Zhang等[9]在實驗和數值上研究了大密度和高黏度的剪切稀化流體中的氣泡動力學問題.其使用Carreau 流變模型作為液相的力學模型,并觀察到氣泡尾部存在一個高黏性區域,這會影響尾隨氣泡的運動,同時剪切稀化流體中流變指數的下降會增加氣泡的上升速度.還有研究表明[10],氣泡變形隨伽利略數(Ga)和厄特沃什數(Eo)的增加或流變指數的降低而增大.氣泡終端速度隨Ga數的增加而增大,隨Eo數和流變指數的增加而減小.Abu-Nab 等[11]研究了表面張力對氣泡半徑、生成時間和氣泡初始速度的影響.胡波等[12]研究了特征時間對氣泡運動的影響,發現特征時間越大,氣泡的終端速度越大.并且當表面張力較小時,氣泡尾部出現黏度盲區.Li 等[13]研究了剪切稀化液中單氣泡的上升,發現氣泡周圍的流線呈現對稱分布,橢球型氣泡的尾部流線更細長.姜韶堃等[14]的研究結果表明,高表觀黏度流體的黏度隨高度呈線性變化,導致氣泡發生規律而可預測的形變,氣泡的寬高比會隨著高度的增加而線性變化,并且氣泡的形狀變化幅度相對較小.Fan 等[15]的研究結果表明,氣泡的瞬時體積和分離體積均隨孔徑減小,隨質量濃度和氣體流量的增加而增大.分離時的長徑比隨孔徑和質量濃度而增大,隨氣體流速的增加而下降.Liu 等[16]研究了非牛頓流體中3 個平行氣泡的相互作用,發現氣泡的水平間距大于臨界間距時會使氣泡發生排斥作用.Wang 等[17]的研究指出,氣泡終端速度容易受壁面效應的影響.還有一些學者對屈服應力和黏彈性流體中氣泡的運動進行了研究.如Xiang 等[18]研究了屈服應力流體中氣泡的形成,發現微通道中液體流量的增加和氣體流量的降低有利于提高氣泡形成的均勻性.Aguirre 等[19]研究了在黏彈性液體中,隨黃原膠濃度的增加,溶液黏度增加而氣泡的速度減小.并發現氣泡周圍產生的渦流是氣泡變形的原因.

除了流變指數、特征時間和表面張力等因素,溫度對非牛頓流體的性質也會產生影響,Rahimi等[20-21]研究了溫度對冪律模型中的稠度指數和流變指數的影響,發現溫度對凝膠的冪律指數的影響不顯著,而稠度指數隨溫度的升高而線性降低,同時在低溫下凝膠的溫度敏感性更強些.Ellahi[22]通過使用兩種與溫度相關的黏度模型,研究了非牛頓納米流體的特性.他們發現黏度參數與溫度呈現強烈的依賴性.

對于處理兩相流問題有許多學者進行了研究并提出了改進方法.Klostermann 等[23]證實使用流體體積方法(VOF)研究氣泡問題具有許多優點,例如能夠保持界面的尖銳性和有界性、保證系統總質量的守恒以及計算快等.Garoosi 等[24]使用了改進版的流體體積方法,對具有不同密度的兩個氣泡的瞬態演變進行了數值研究.Reza 等[25]使用分段線性界面計算方法(PLIC)跟蹤了兩相流界面,通過連續表面力(CSF)模型計算表面張力,并在數值模擬中取得了兩個連續斜坡處氣泡運動的最大速度.Nahed 等[26]基于PLIC-VOF 方法,提出了一種新的曲率估計方法,大大降低了偽流動的影響,使得曲率計算精度更高.

目前對于溫度敏感性凝膠的研究還主要集中在本構方程構建和物性實驗測量等方面,對于氣泡動力學研究還相對較少.本工作采用流體體積方法研究了在靜止凝膠流體中,溫度、特征時間、流變指數和表面張力對氣泡形變和運動的影響,重點關注氣泡的縱橫比、速度、重心位置的變化趨勢,以期為凝膠推進劑加注過程中的除氣工藝提供研究基礎.

1 物理模型和計算工況

1.1 物理模型

本文研究的氣泡運動速度較小,三維效應不顯著,因此為了節省運算時間,采用二維模型進行仿真研究.如圖1 所示,計算域為一個二維矩形區域.Pang 等[10]采用二維模型仿真研究了氣泡在非牛頓流體中的上升運動,與實驗結果誤差為12%.值得指出的是本文所簡化的模型為二維無限大平板,與三維球形氣泡存在差異.氣泡初始位置位于A點受浮力驅動向上運動,其中L足夠長,可以確保氣泡的運動得到充分發展.由文獻[27]可知計算域寬度的一半與氣泡直徑的比值超過10 時,氣泡的終端速度趨于穩定,此時可忽略壁面效應的影響.本文研究的氣泡半徑較小,W/(2d)=12.5,所以壁面對氣泡的運動影響可以忽略不計.矩形計算區域的寬度為W=0.15 m,高度為L=0.24 m,初始氣泡中心A點與底部壁面距離為h=0.04 m,定義氣泡初始位置為原點.網格尺度為0.000 2 m,一個氣泡內約有700 個網格.

圖1 示意圖(非按比例)Fig.1 Schematic (not to scale)

假設氣泡的初始形狀為圓形,直徑0.006 m,氣泡的上升運動主要受重力、黏性力和表面張力的影響,其中重力加速度的方向向下.這些作用力對于氣泡運動的影響可用伽利略數(Ga) 和厄特沃什數(Eo)進行表征.其中,Ga數表示重力與黏性力的比值,Eo數表示重力與表面張力的比值.其表達式為

其中,ρl表示液相密度,ρg表示氣相密度,g表示重力加速度,d表示氣泡初始直徑,η0表示液相零剪切黏度,σ 表示表面張力,θ 表示溫度.

1.2 邊界條件

在模擬中,上端邊界為壓力出口邊界條件(pressure outlet),而下端邊界和左右側面邊界為無滑移壁面邊界條件(wall).

2 控制方程和數值方法

2.1 控制方程

本文假設氣液兩相的流動均為不可壓縮的且為層流,并忽略由溫度變化引起的流體密度變化.Navier-Stokes 方程的守恒形式[28],包括連續性方程和動量守恒方程的公式如下

其中,U為速度,ρ 為流體局部平均密度,? 表示張量積,p為壓力,τ 為應力張量,F為運用連續表面張力模型得到的體積力源項,? 為哈密頓算子.其中,應力張量的公式如下

其中,μ 為流體當地平均黏度,D為應變率張量,其公式如下

流體局部平均密度 ρ 和平均黏度 μ 的計算如下

其中,α 為整個計算域的相函數,其定義見2.3 節;下標l 和g 分別表示液相和氣相.利用連續表面張力模型[29],我們可以將表面張力當作體積力處理,計算體積力的方法如下

式中,σ 為表面張力系數,κ 為界面曲率,計算公式如下

本文研究的凝膠溶液與氣體的黏度比值較大,約為 1 05,而凝膠溶液與氣體的密度比為1 03.研究指出如果氣液之間的黏度和密度比太大,則計算收斂非常困難[30].因而對于壓力方程采用共軛梯度求解器和多重網格預條件器進行求解,為了保證數值求解的收斂性和穩定性,控制時間步長,使得庫朗數最大在0.05 左右.同時,采用PIMPLE 算法來求解速度場和壓力場的耦合問題,外部迭代50 次,內部迭代3 次.在網格上使用分段線性界面重建來處理流體界面上的不連續性,以實現數值求解的連續性.

2.2 液相本構方程

凝膠的剪切稀化特性采用Carreau-Yasuda 模型描述,根據局部剪切速率,可得到凝膠的當地表觀黏度.Carreau-Yasuda 模型如下

重點考慮了溫度對凝膠推進劑黏度的影響,并根據Andrade-Eyring 法則[31],使用Reynolds 模型來表示溫度影響項[32],溫度敏感剪切變稀凝膠的黏度本構方程如下

式中,η∞和 η0分別為無窮大剪切率和零剪切率時的黏度,λ為時間常數,b控制了黏度-剪切率曲線初始階段的形狀,n為流變指數,c1和c2為溫度指數.γ˙ 為剪切速率,計算如下

剪切速率和溫度是影響黏度的兩個獨立參數.由課題組實驗測量與參數擬合獲得[33],具體的本構方程參數為 η0=71.217 3 Pa·s,η∞=0.027 1 Pa·s,λ=0.158 2,b=1.381 4,n=0.396 4,c1=-0.01 1,c2=0.061 1.

2.3 流體體積法

兩相流問題中,界面的復雜演化是宏觀尺度兩相流問題的一個主要特征[34].采用VOF 方法來跟蹤流體表面.該方法利用流體體積與網格體積的比值來計算VOF 函數的值,該函數通過界面壓縮方法保持尖銳的界面,并且在沒有重新初始化步驟的情況下保持質量守恒.通過求解VOF 函數所滿足的運輸方程,可以獲得全流場的流體體積函數分布情況,并進而構建運動界面.通過追蹤界面,可以獲得其附近所有的流場信息.流體運動時,氣-液的交界面也變化,所以需要更新每個單元對應的體積函數,相函數的運輸方程[35]是VOF 的核心,表達式為

其中,u為速度矢量;α 為整個計算域的相函數,其值在氣泡內部為0,在氣泡外部為1,在界面處為0～1 之間.下式用于計算氣相的局部體積分數

其中,下標s 是網格編號;As是編號為s 的網格單元的面積;ψs是界面追蹤的直觀指標,同時也是界面構造的關鍵參數.

此外,還采用PLIC 方法對界面進行重構,通過在單個網格中劃分斜線線段來近似兩相界面,斜線的斜率通過附近的相分數來確定,具有誤差小、精度高的優點.通過求解單個網格內氣泡的局部體積分數梯度,可以獲得界面的法向量.計算界面法向量的方法如下

在每個網格中,下一時刻氣泡的局部體積分數可以通過計算當前時刻網格內氣泡的局部體積分數與這段時間流入網格內的氣泡體積分數之和來獲得.通過更新下一時刻網格內的氣泡局部體積分數,就可以對整個流場中的氣泡界面進行重構.

3 結果分析和討論

3.1 不同溫度下的氣泡運動特性

氣泡的運動速度和周圍液體的流動狀態會影響氣液兩相之間的傳質傳熱效率[36].氣泡的變形程度受到黏性力和表面張力的影響,黏性力或表面張力越大,氣泡變形程度越小.研究將體積分數為0.5 的區域定義為氣泡區域.圖2(a)顯示了不同溫度下氣泡穩定后的輪廓,由于凝膠液體黏度較大,所以本文中所有氣泡都沒有明顯的變形,氣泡形狀更接近球形.圖2(b)和圖2(c)顯示了不同時刻下氣泡的運動軌跡和重心位置,可以看到溫度升高會使氣泡運動得更遠,并在豎直方向基本保持直線運動.氣泡的變形是由于氣泡的上升過程中,其頂端與底端壓力的分布不均導致的.如圖2(d)所示,溫度升高會增加氣泡的變形.由于氣泡的形狀均表現為橢圓形,因此氣泡的變形程度主要表現在縱橫比上.縱橫比是衡量氣泡變形程度的重要參數,其值為氣泡最大橫向距離與最大縱向距離之比[37].

圖2 不同溫度下氣泡的運動情況Fig.2 The movement of bubbles at different temperatures

氣泡形狀的變化會對氣泡的動力學特性產生重要的影響[38],尤其是氣泡的尾渦會隨氣泡形狀的變化而不斷變化,并最終達到穩定狀態.圖3 顯示了氣泡中心參考系下的流線圖,展示了氣泡的流動情況.氣泡的變形會導致流線的不同,可以看到在氣泡內部出現兩個環流區,而氣泡的尾流區沒有任何額外的環流區.由于內部環流的慣性作用,環流會被底部表面切為兩個部分,即內部環流和外部渦流.由于變形氣泡的曲率較小,氣泡尾部沒有形成尾渦[39].圖4展示了不同溫度下氣泡上升速度的變化情況.氣泡上升的初始階段浮力占主導作用,導致氣泡的縱向速度急劇增加[40].隨著氣泡的不斷上升,在高黏度的流體中氣泡所受到的阻力逐漸增大,使得氣泡上升的加速度變小,直到氣泡達到一個穩定狀態.此時氣泡的速度即為其終端速度.隨著溫度升高,液相的表觀黏度降低,氣泡周圍液體的阻力也隨之降低,導致氣泡的終端速度隨著溫度升高而增加.圖4(b)顯示了氣泡上升速度與縱橫比的關系.可以看到氣泡的上升速度與縱橫比呈正相關.氣泡速度[41]的計算公式為

圖3 不同溫度下氣泡周圍的流線圖Fig.3 Streamline diagram around bubbles at different temperatures

圖4 不同溫度下氣泡的速度與縱橫比Fig.4 Velocity and aspect ratio of bubbles at different temperatures

式中,Ω 為氣泡區域的面積,u為每個網格上讀取的速度.

3.2 不同溫度下氣泡周圍液相黏度及剪切速率分布

由于凝膠流體的黏度較高,氣泡在其中運動時軌跡和形狀相對穩定,從而容易形成穩定的剪切流場.不同氣泡的終端速度差異主要是由于它們受到的阻力不同[42].摩擦阻力與液相的黏度相關,圖5 上方顯示了氣泡周圍液相的表觀黏度分布情況.為了方便比較,將液相的表觀黏度進行了最大值歸一化處理,4 條等高線分別為0.98,0.95,0.9 和0.8.剪切速率的增加會降低剪切變稀流體的黏度,因此當氣泡上升時,氣泡周圍的液體黏度會降低,即處于高剪切區.同時,氣泡周圍液相低表觀黏度區的增大會減小氣泡上浮過程所受到的摩擦阻力,從而提高氣泡的終端速度.

圖5 不同溫度下液體中的表觀黏度(上)和剪切速率(下)Fig.5 Apparent viscosity (top) and shear rate (bottom) in liquids at different temperatures

圖5 下方為剪切速率分布的云圖,根據圖中的分布情況可以看出,在剪切速率較低的區域,液相表觀黏度較高,反之則較低.此外,在氣泡的赤道處出現了一條腰線,氣泡頂部和底部均出現了高剪切速率區,并從氣泡上下兩部分向外輻射.在腰線處液體的表觀黏度最高,剪切速率最低.

3.3 流體性質對氣泡運動的影響

當流體的流變指數n=1 時,該流體為牛頓流體,其黏度與剪切速率無關.當n＜ 1 時,流體表現為剪切變稀流體.流變指數n控制黏度隨剪切速率快速變化期間的斜率.圖6 展示了n=0.2,n=0.4,n=0.6 和n=0.8 這4 種情況下氣泡周圍液相的表觀黏度和剪切速率分布.結果表明,當n值較小時,液相的剪切稀化效應更明顯,說明流體更容易流動,氣泡周圍液體的表觀黏度降低,剪切速率增加.

圖6 不同流變指數下的表觀黏度(左)和剪切速率(右)Fig.6 Apparent viscosity (left) and shear rate (right) at different rheological indices

圖7(a)顯示了流變指數n對靠近氣泡底部的液體黏度的影響.隨著流變指數的減小,氣泡底部附近的液體黏度下降更多.同時,剪切稀化效應會導致液體的黏度降低,動能耗散減小,從而也會使氣泡速度增加.圖中還顯示了氣泡在牛頓流體n=1 中上升的結果.如圖7(b)說明了流變指數n對靠近氣泡下方的液體速度的影響.越靠近氣泡下方液體速度越大,其中由于剪切稀化效應,其值隨著流動指數的減小而增加.

圖7 流變指數對氣泡底部的黏度及速度的影響Fig.7 Effect of rheological index on viscosity and velocity at bubble bottom

當氣泡向上運動時,液體被推向氣泡的下方,形成一個向下的流動.這個流動使得氣泡周圍的液體受到強烈的剪切應力,從而引起液體的黏度降低,同時也使液體的流動速度增加.另外,浮力的作用使得氣泡周圍的液體向外擴散,這進一步降低了液體的黏度.

圖8(a)和圖8(b)分別顯示了在剪切稀化流體中,不同流變指數下氣泡縱橫比和重心隨時間的變化.從圖中可以觀察到,在剪切稀化效應越明顯的流體中,氣泡的縱橫比變化越小.液體中低表觀黏度區的增大是由于剪切稀化效應導致的,這使得氣泡受到的黏性阻力減小,其上升速度變快,因此氣泡的重心位置也逐漸升高.

圖8 流變指數對氣泡的縱橫比及重心位置的影響Fig.8 Effect of rheological index on the aspect ratio and the position of the center of gravity of bubbles

氣泡終端速度會受到液相流體性質的影響,如圖9(a)所示,隨著時間的推移,氣泡速度快速增加并逐漸穩定.在剪切稀化效應較強的流體中,由于流體的剪切作用,氣泡周圍液相的黏度會大幅下降,從而減小了氣泡受到的黏性阻力,因此在剪切稀化液體中氣泡更容易上升.氣泡的終端速度也越大.

圖9 流變指數及特征時間對氣泡速度的影響Fig.9 Influence of rheological index and characteristic time on bubble velocity

同時,特征時間λ是描述剪切稀化液體特性的重要參數.當液相受到擾動時,液相中的高分子鏈受到擾動,平衡態會遭到破壞,特征時間λ是描述液相從非平衡態恢復到平衡態所需要的時間,當λ=0 時,流體為牛頓流體.在本研究中,選擇了λ=0.05 s,0.158 2 s,0.25 s,0.5 s,1 s 這5 種情況.從圖9(b)可以看出,氣泡的終端速度隨著λ的增加逐漸增加,在t=0.4 s 之前,氣泡終端速度增幅較大,而在t=0.4 s 后,氣泡終端速度增幅較小.此外,當特征時間較大時,其對氣泡速度的影響也更為明顯.

圖10 顯示了氣泡的運動軌跡隨著特征時間的增加而變化,氣泡在縱向拉長.這是因為特征時間越大,λ越大,液相受到擾動后恢復到平衡態所需的時間越長,剪切稀化效應更強,氣泡形狀越不穩定,氣泡被拉伸和扭曲的程度更高.

圖10 不同λ 下的氣泡運動軌跡Fig.10 Bubble trajectory under different λ

根據圖11(a)和圖11(b),可以看到在不同的溫度和流變指數下,氣泡終端速度和縱橫比均有相應的變化.隨著液體流變指數的增加,氣泡周圍的阻力會增加,導致氣泡終端速度下降,同時氣泡的形變也會增加.此外,如前所述,溫度對液體的黏度也會產生影響,從而間接地影響氣泡周圍的摩擦和阻力大小.因此,隨著溫度升高,氣泡的終端速度也會增大,并且氣泡的變形程度也更為明顯.在低溫下流變指數對氣泡縱橫比的影響要更小些.

圖11 3 種溫度下氣泡的縱橫比和速度隨流變指數與特征時間的變化情況Fig.11 Variations of aspect ratio and velocity of bubbles with rheological index and characteristic time at three temperatures

根據圖11(c)所示,λ的增加加速了液體從牛頓流體到剪切稀化流體的轉變.當溫度保持不變時,液體的λ減小會使其黏性下降,氣泡周圍的液體阻力隨之下降,使得氣泡的縱橫比更大.同時隨溫度的升高,氣泡的縱橫比在均勻增加.由圖11(d)知氣泡的終端速度隨著λ 的增加逐漸增加,并且高溫能夠促進氣泡的上升,使得氣泡的終端速度增加.隨λ的增加,溫度對氣泡終端速度的影響逐漸變大,氣泡速度增加的幅度也隨之變大.

在大的Eo數條件下,由于表面張力相對較弱,氣泡會具有較大的變形.如圖12 所示,在本文研究的Eo=15,100,950 這3 種情況下,隨著計算時間的增加,氣泡底部繼續向上凹陷并向外水平擴展.當表面張力較大時,氣泡可以抵抗射流的影響并保持形狀穩定,使得射流不會突破氣泡邊緣,射流在底部界面沿徑向方向流動[43].而當表面張力較小時,在氣泡內外壓力差的作用下,氣泡的底部向上凹陷,但并不會在氣泡外部形成渦流.

圖12 不同Eo 數下的流線圖Fig.12 Streamline diagrams under different Eo number

如圖13 所示,隨著表面張力的減小,氣泡底部出現黏度盲區,這會使氣泡的上升運動受到一定程度的阻礙[44].從圖14 可以看到,當表面張力減小時,氣泡的終端速度會增加,但當Eo數超過100 時,這時再減小表面張力也不會使氣泡的終端速度發生明顯的變化.

圖13 不同Eo 數下液體的表觀黏度Fig.13 Apparent viscosity of liquid at different Eo number

圖14 不同Eo 數下氣泡的上升速度Fig.14 Rising speed of bubbles at different Eo number

4 結論

本文探究了凝膠流體中氣泡的動力特性.通過采用VOF 數值模擬方法,研究了凝膠流體中不同溫度、流變指數、特征時間和表面張力等參數對氣泡動力特性的影響,并得出如下結論.

(1)隨著溫度升高,凝膠的表觀黏度降低,剪切速率增加,氣泡所受阻力減小,導致氣泡終端速度增加,且變形程度增加.

(2)改變凝膠流體的流變特性會對氣泡的縱橫比、重心位置和終端速度產生不同程度的影響.液體剪切稀化趨勢的增強使氣泡周圍的液體黏度顯著降低,進而使氣泡在剪切變稀液體中更容易上升.特征時間對氣泡速度的影響比較顯著,特征時間的增大會增加氣泡速度,并且增加的幅度比改變流變指數造成的影響要大.

(3)溫度、流變特性和特征時間的變化對氣泡變形影響相對較小,但表面張力對氣泡底部的變形影響較大.特別是在高Eo數下,即低表面張力下,氣泡在上升過程中受到的阻力較小,因此氣泡底部會向上凹陷.

綜上所述,研究結果有助于深入理解氣泡在不同流體特性下的運動和變形特性,能為未來凝膠使用過程的優化設計提供重要依據.