二維平行平板內(nèi)分散氣泡對液相表觀黏度影響的數(shù)值研究

徐磊，龐明軍

（常州大學機械工程學院，江蘇 常州 213164）

二維平行平板內(nèi)分散氣泡對液相表觀黏度影響的數(shù)值研究

徐磊，龐明軍

（常州大學機械工程學院，江蘇 常州 213164）

泡狀液體廣泛存在于各種工業(yè)領(lǐng)域，因此深入了解泡狀液體的物理特性對產(chǎn)品質(zhì)量改進和過程強化具有重要意義。本文為了理解氣泡對液相表觀黏度的影響，揭示液相表觀黏度變化的原因，利用流體體積函數(shù)（VOF）與動網(wǎng)格結(jié)合的方法對二維平行板模型內(nèi)氣泡對液相表觀黏度的影響進行了深入的數(shù)值研究；分別研究了相同體積分數(shù)、不同毛細數(shù)工況下，氣泡加入對液相表觀黏度的影響。結(jié)果表明：毛細數(shù)和氣泡體積分數(shù)對液相的表觀黏度具有巨大的影響作用；相同的體積分數(shù)下，當毛細數(shù)較大（Ca＞1）時，氣泡的加入會使液相的相對黏度減小，體積分數(shù)越大，液相的相對黏度減小得越嚴重；另外，相同的體積分數(shù)下，當毛細數(shù)較小（Ca＜1）時，氣泡的加入會使液相的相對黏度增大，體積分數(shù)越高，液相的相對黏度增大得越明顯。

泡狀液體；分散氣泡；相對黏度；流體體積函數(shù)；動網(wǎng)格；二維平行平板

泡狀液體廣泛存在于化工、食品、流動輸送、生物醫(yī)學等工業(yè)領(lǐng)域。已有的研究表明，氣泡的加入會極大地影響工業(yè)過程的質(zhì)量、熱量和動量傳遞效率，甚至會對泡狀設施的運行效率和產(chǎn)品質(zhì)量帶來巨大的影響［1-11］。為了有效地設計泡狀設施和高效地控制泡狀設施的運行過程，深入理解泡狀液體的水動力學特性是非常必要的。特別是理解氣泡加入對液相流變特性的影響，對于理解氣泡湍流減阻機理和改進食品工業(yè)的產(chǎn)品質(zhì)量具有重要的意義。

目前，國外研究者針對泡狀液體的流變特性已開展了一定的實驗和數(shù)值研究。RUST和 MICHAEL［1］基于旋轉(zhuǎn)流變儀原理測量了氣泡對玉米糖漿黏度的影響。YUICHI和HIROSHI［2］利用落球黏度計原理詳細測試了氣泡瞬態(tài)變形對硅油黏度的影響。LLEWELLIN等［3］學者在實驗研究的基礎上，通過理論分析給出一組泡狀流流變特性的本構(gòu)模型。文獻［4-8］報道了氣泡對巖漿流變特性影響的研究。文獻［9-10］給出了低毛細數(shù)下，泡狀流相對黏度的本構(gòu)方程。文獻［7，11］運用數(shù)值方法對泡狀流的流變特性開展了研究。

綜上所述，發(fā)現(xiàn)盡管針對泡狀液體的流變特性開展了一定的研究，但限于問題的復雜性，仍未獲得 一個很好的機理解釋，不同的學者給出不同的結(jié)果。為了進一步理解氣泡存在時是如何影響液相的表 觀黏度，并希望獲得一些液相表觀黏度變化的原因，本文采用二維平行板模型，對分散氣泡對液相表觀黏度的影響進行了數(shù)值模擬。

1 幾何模型和計算條件

1.1 幾何模型、物性參數(shù)與邊界條件

鑒于二維平行平板流動（即庫埃流動模型）簡單，在低雷諾數(shù)下，能給出理論解，且結(jié)果便于分析和對比，為此選用二維平行平板作為計算的幾何模型，如圖1所示。其中平板長度為L=1m，板間距為h=0.5m；計算時，將氣泡置于二維平行平板中心，氣泡半徑為r=0.1m、0.13m、0.15m。為了減小計算工作量，減小計算區(qū)域，在左右邊界設置為周期邊界條件；讓下平板靜止不動，而上平板以恒定的速度u=0.1m/s向右作勻速運動。為了便于計算收斂，將氣液相的密度均設置為ρ=0.1kg/m3，液相的動力黏度設置為μl=1000Pa·s，氣相的動力黏度設置為μb=100Pa·s。

1.2 網(wǎng)格劃分

對于流體體積函數(shù)（VOF）方法而言，選用非結(jié)構(gòu)化網(wǎng)格，計算精度較高、便于界面重構(gòu)；且非結(jié)構(gòu)化網(wǎng)格由于網(wǎng)格單元形狀的多樣性以及單元之間連接關(guān)系的自由性、適應強，為此本文選用三角形網(wǎng)格進行網(wǎng)格劃分。其中，網(wǎng)格節(jié)點數(shù)為2145，網(wǎng)格單元數(shù)有4096，網(wǎng)格具體劃分情況見圖2。

圖1 計算幾何模型

圖2 網(wǎng)格劃分

1.3 計算工況

目前的研究表明，氣泡影響液相表觀黏度的因素可概括為毛細數(shù)（Ca）和體積分數(shù)（Φ）兩個量綱歸一化參數(shù)。為此，根據(jù)目前的計算模型和計算能力，依據(jù)這兩個參數(shù)設計了如下計算工況，詳見表 1。體積分數(shù)的改變是通過調(diào)整氣泡半徑的大小實現(xiàn)的，毛細數(shù)的大小，本文是通過改變表面張力系數(shù)實現(xiàn)的。計算時，為了便于計算收斂（即便于界面重構(gòu)），氣泡半徑選取得較大。

其中毛細數(shù)是衡量氣泡變形大小的一個參數(shù)［1-5］，具體指的是氣泡所受剪切力與表面張力的相對大小。計算表達式詳見式（1）。

式中，μ0為純液體的動力黏度，Pa·s；為剪切率，s-1；r為沒有變形時氣泡的半徑，m；表示的是液相的表面張力系數(shù)，N/m。

表1 計算工況（含對應工況下的上平板所受剪應力）

2 控制方程和計算方法

2.1 控制方程

因目前氣液兩相的運動速度較小，所以計算時均不考慮其可壓縮性。氣液兩相混合流體的控制方程表述如式（2）、式（3）。

式中，ρ是氣液混合相的密度，kg/m3；p為壓強，Pa；t是時間，s；u為速度矢量，m/s；μ是黏度，Pa·s；F表示的是表面張力，N。

2.2 VOF方法

精確描述氣液兩相的運動界面是研究氣泡行為的關(guān)鍵。由于 VOF 方法追蹤的是網(wǎng)格中的流體體積，不是追蹤流體質(zhì)點的運動，因而具有計算速度較快和界面追蹤精度高等優(yōu)點［12-16］。因此，本文首先考慮采用VOF方法追蹤氣泡-流體界面。VOF方法的核心是相函數(shù)的輸運方程，所有的界面構(gòu)造和流體輸運都基于此方法進行計算，其表達式為

式（4）。

式中，u為速度，m/s；f為定義在整個計算區(qū)域內(nèi)的相函數(shù)。

公式（4）可以看作f的守恒輸運公式，以二維流場為例，此時若將流場離散為網(wǎng)格單元，將f值在單個網(wǎng)格面積上求積分，該積分值與當前網(wǎng)格單元面積之比即為流體的局部體積分數(shù)（ψ），即式（6）。

式中，i是網(wǎng)格單元編號；Ai是編號為i的網(wǎng)格單元面積。該公式描述了網(wǎng)格局部體積分數(shù)ψ與相函數(shù)f之間的關(guān)系，其中f是流體計算的對象，而ψ是界面追蹤的直觀指標以及界面構(gòu)造方案的關(guān)鍵。

當計算單元中是液體和氣體時，方程中的流體物性為相應液體或氣體的物性。當計算單元內(nèi)包含兩相界面時，流體物性按照兩相體積分數(shù)的加權(quán)平均進行計算，即

2.3 表面張力

氣泡對液相的影響是通過表面張力實現(xiàn)的，要求解動量方法，必須計算出動量方程（3）中的表面張力項。

本文采用連續(xù)表面張力模型（即CSF模型）來計算作用在氣液兩相界面的上力，其計算表達式為式（8）。

式中，k為流體自由面的局部曲率；δ為狄拉克分布函數(shù)；n為自由面的局部單位法向量。

2.4 界面重構(gòu)

界面重構(gòu)是VOF模型最關(guān)鍵的部分，這里應用分段線性界面計算方法（PLIC）［13］，對氣液兩相界面進行重構(gòu)。該界面重構(gòu)方法的特點就是在單個網(wǎng)格內(nèi)用直線段近似氣液兩相的界面。首先確定界面的法向量n。界面的法向?qū)嶋H上就是相函數(shù)f的梯度方向，可以通過求解網(wǎng)格單元中氣泡局部體積分數(shù)（ψ）的梯度獲得［15］。界面法向單位法向量的計算公式如式（9）。

然后根據(jù)法向和當前網(wǎng)格中的氣泡相的局部體積分數(shù)計算出當前氣液界面線的確切位置。最后將所有網(wǎng)格的氣泡相的局部體積分數(shù)加上流入的流體再減去流出的流體就是每個網(wǎng)格中下一時刻新的氣泡相的局部體積分數(shù)，依據(jù)新的氣泡相的局部體積分數(shù)就可以重新對整個流場的自由面進行界面重構(gòu)。

2.5 數(shù)值求解過程

對于目前的計算，因上平板作勻速運動，所以采用動網(wǎng)格進模型進行計算。選用層流模型進行計算，求解時速度與壓力場的耦合采用SIMPLE算法；單元中心的變量梯度采用Least-Squares cell-Based，壓力的離散采用 PESTO格式，動量方程采用二階迎風格式進行離散，體積分數(shù)采用Geo-Reconstruct格式進行離散；時間步長Δt=0.001，時間的推進采用一階隱式格式。

3 數(shù)值模擬結(jié)果分析

為了驗證計算結(jié)果的準確性，先對單相流動進行了模擬，將模擬結(jié)果與理論結(jié)果進行了對比，詳見圖 4?？梢钥闯觯壳暗挠嬎憬Y(jié)果與理論結(jié)果吻合的很好，表明計算結(jié)果是可信的。為了觀察計算結(jié)果是否取得統(tǒng)計穩(wěn)態(tài)，計算兩相流動時，考察了上平板上的剪切力隨計算時間的變化情況，如圖5所示?？梢钥闯觯?000次以后，上平板上的剪切力基本趨于穩(wěn)定，表明計算結(jié)果可以進行統(tǒng)計處理。本文選取迭代5000次以后的值進行平均，從而得到上平板的剪切力。然后再利用牛頓流體剪應力的計算公式反推得到液相的表觀黏度。

圖4 單相流的速度計算曲線

圖5 監(jiān)控體積分數(shù)為Φ=0.063與毛細數(shù)Ca=10的上平板上的剪切力

3.1 氣泡加入對液相表觀黏度的影響

為了研究氣泡加入對液相表觀黏度的影響，圖6給出不同工況下液相表觀黏度的變化情況。其中橫坐標為毛細數(shù)，縱坐標為液相的相對黏度（為加入氣泡后液相的表觀黏度除以純液體對應狀態(tài)下的黏度）。為了對比，圖6也給出了文獻［8］中的實驗結(jié)果，可以看出目前的計算結(jié)果與文獻［8］的實驗結(jié)果定性一致，定量的差別或許是因目前計算時給出的假設引起的。

圖6 相對黏度變化情況

從圖6可以看出，當毛細數(shù)小于0.5時（Ca＜0.5），隨著毛細數(shù)的減小，氣泡加入使液相的相對黏度增大，當毛細數(shù)小到一定程度時（對于目前的計算Ca≈0.2），液相相對黏度趨于穩(wěn)定、不再隨著毛細數(shù)的減小而再增大；然而，當毛細數(shù)大于 0.5時（Ca＞0.5），隨著毛細數(shù)的增大，氣泡加入使液相的相對黏度減小，當毛細數(shù)大到一定程度時（對于目前的計算Ca≈10），液相相對黏度趨于穩(wěn)定、不再隨著毛細數(shù)的增大而再減小。而且在相同的毛細數(shù)下，氣泡的體積分數(shù)越大，氣泡對液相相對黏度的影響就越大。

由此看來，氣泡對液相表觀黏度的影響與氣泡的體積分數(shù)和毛細數(shù)有直接的關(guān)系。而毛細數(shù)的大小又直接決定了氣泡形變的大小，所以氣泡形狀的改變又會對液相相對黏度的變化產(chǎn)生直接的影響。為此，選取了最大體積分數(shù)（Φ=0.14）下氣泡形狀隨毛細數(shù)的變化情況，詳見圖7。圖7為計算時間為5s時的氣液相相分布云圖，其中藍色區(qū)域表示的是氣泡，紅色為液相所占區(qū)域?？梢钥闯觯S著毛細數(shù)的增大，氣泡由圓形逐漸變?yōu)榧氶L的橢圓形，且橢圓形氣泡的長軸與平板法向的夾角與最大剪切率與法向的夾角基本一致。

當毛細數(shù)較小時（對于目前的計算 Ca＜0.2），氣泡所受的剪切力小于其所受的表面張力，氣泡保持為圓形，這時液相的相對黏度由于氣泡的加入而增大。當毛細數(shù)進一步減小時，氣泡的形狀不會發(fā)生變化、一直保持為球形，所以對應的液相的相對黏度也不會再增大、而趨于一個定值。然而，當毛細數(shù)較大時（對于目前的計算Ca＞10），氣泡所受的剪切力大于表面張力，由于強的剪切作用，氣泡不能再保持球形而變?yōu)榧氶L的橢圓形，這時液相的表觀黏度由于氣泡的加入而減小。因這時氣泡長軸與平板法向的夾角已近似等于最大剪切率與法向的夾角，當毛細數(shù)再增大，橢圓形氣泡長軸與法向的夾角也不會變化，所以液相的相對黏度也趨于一個定值、不會隨著毛細數(shù)的增大而減小。

圖7 氣泡體積分數(shù)Φ=0.14在5s時氣液相分布云圖

3.2 氣泡加入對液相流場的影響

為了理解毛細數(shù)（即氣泡形狀）對液相表觀黏度的影響，下面分析了氣泡加入對液相流場的影響情況。圖8給出最大體積分數(shù)為Φ=0.14下毛細管為0.1和1時的壓力場分布，圖9給出速度矢量場和渦量場的分布情況。可以看出，當毛細數(shù)為 0.1時，由于較大的表面張力，氣泡幾乎保持為球形，氣泡周圍的壓力小于氣泡內(nèi)部的壓力，在氣泡周圍最大剪切方向上產(chǎn)生大面積的負壓區(qū)，對應的在負壓區(qū)域出現(xiàn)大面積的渦量區(qū)（如圖9所示），導致氣泡與液體在此區(qū)域的黏性耗散增加，所以宏觀上表現(xiàn)為液相相對黏度的增加。

圖8 氣泡體積分數(shù)Φ=0.14的壓力場

當毛細數(shù)為1時，由于剪切力大于表面張力，氣泡發(fā)生嚴重的變形，在流體流動方向［如圖 8（b）中箭頭所示］上，氣泡的左右兩側(cè)出現(xiàn)了小面積的負壓區(qū)，對應的在此區(qū)域也出現(xiàn)了一定的渦量區(qū)；負壓區(qū)沒有出現(xiàn)在最大剪切方向上，所以在該方向上沒有渦運動（如圖9所示），即氣泡形變提供了自由滑移面，所以氣泡加入使黏性耗散增加的不明顯。反而，由于氣泡形變后沿著最大剪切方向呈現(xiàn)為細長橢圓形，而氣體的黏度又小于液體的黏度，即低黏度流體占據(jù)了大面積的高剪切率區(qū)域，所以形變的氣泡導致黏性耗散大幅度減小?？傊?，對于形變的氣泡而言，其黏性耗散減小的程度大于黏性耗散增加的程度，最終宏觀表現(xiàn)為液相的相對黏度減小。

為了直觀得觀察不同毛細數(shù)下，氣泡對液相流場的影響情況，圖10給出了不同毛細數(shù)下流場流線分布情況?？梢钥闯?，與單液體相相比，毛細數(shù)較小時［如圖10（b）所示］，氣泡形變較?。ㄒ簿褪钦f氣泡幾乎保持為圓形），氣泡加入，使流線彎曲的越厲害，導致黏性耗散增加，表現(xiàn)為液相相對黏度增大；反之，毛細數(shù)較大時［如圖10（d）所示］，氣泡形變較大、為流體的流動提供了滑移面，氣泡加入使流線的彎曲不明顯，導致總體黏性耗散減小，表現(xiàn)為液相相對黏度減小。這與文獻［1，5］給出的解釋相符。

圖9 氣泡體積分數(shù)Φ=0.14的速度場和渦量場

圖10 不同毛細數(shù)下流線分布

4 結(jié) 論

本文通過對二維平行板模型內(nèi)分散氣泡對液相表觀黏度的影響進行了詳細的數(shù)值模擬，可以得到如下幾點結(jié)論。

（1）當毛細數(shù)小于0.5時，氣泡所受的剪切力小于表面張力，氣泡幾乎保持為圓形，氣泡的加入使液相的相對黏度增大。當毛細數(shù)進一步減小時，由于氣泡的形狀不會發(fā)生變化、一直保持為球形，所以對應液相的相對黏度也不會再增大，而趨于一個定值。而且在相同的毛細數(shù)下，氣泡的體積分數(shù)越大，氣泡對液相相對黏度增加的越大。

（2）當毛細數(shù)大于0.5，氣泡所受的剪切力大于表面張力，在剪切作用下，氣泡變?yōu)榧氶L的橢圓形，氣泡的加入導致液相表觀黏度的減小。這時氣泡長軸與平板法向的夾角已近似等于最大剪切率與法向的夾角，當毛細數(shù)再增大，橢圓形氣泡長軸與法向的夾角也不會變化，所以液相的相對黏度也趨于一個定值，不會隨著毛細數(shù)的增大而減小。然而隨著體積分數(shù)的增加，液相表觀黏度的減小幅度增大。

（3）在低毛細數(shù)下，由于較大的表面張力，氣泡幾乎保持為球形，氣泡周圍的壓力小于氣泡內(nèi)部的壓力，在氣泡周圍最大剪切方向上產(chǎn)生大面積的負壓區(qū)，對應的在負壓區(qū)域出現(xiàn)大面積的渦量區(qū)，并且氣泡的加入使流線彎曲嚴重，導致黏性耗散增加，所以表現(xiàn)為液相相對黏度增大。

（4）在高毛細數(shù)下，由于氣泡所受的剪切力大于表面張力，氣泡發(fā)生嚴重的變形，在流體流動方向上，氣泡的左右兩側(cè)出現(xiàn)了小面積的負壓區(qū)，對應的在此區(qū)域也出現(xiàn)了一定的渦量區(qū)；負壓區(qū)沒有出現(xiàn)在最大剪切方向上，所以在該方向上沒有渦運動，即氣泡形變提供了自由滑移面。所以形變的氣泡導致黏性耗散大幅度減小。表現(xiàn)為液相相對黏度減小。

［1］ RUST A C，MICHAEL M. Effects of bubble deformation on the viscosity of dilute suspensions［J］. J. Non-Newtonian Fluid Mech.， 2002，104：53-63.

［2］ YUICHI M，HIROSHI O. Increase of effective viscosity in bubbly liquids from transient bubble deformation［J］. Fluid Dynamics Research，2008，40：565-575.

［3］ LLEWELLIN E W，MADER H M，WILSON S D R. The rheology of a bubbly liquid［J］. Proc. R. Soc. A，2002，458：987-1016.

［4］ MANGA M，LOEWENBERG M. Viscosity of magmas containing highly deformable bubbles［J］. Journal of Volcanology and Geothermal Research，2001，105：19-24.

［5］ MADER H M，LLEWELLIN E W，MUELLER S P. The rheology of two-phase magmas：a review and analysis［J］. Journal of Volcanology and Geothermal Research，2013，257：135-158.

［6］ LEJEUNE A M，BOTTINGA Y，TRULL T W，et al. Rheology of bubble-baring magmas：reply to a critical comment by Spera and Stein［J］. Earth and Planetary Science Letters，2000，175：333-334.

［7］ MANGA M，CASTRO J，CASHMAN K V，et al. Rheology of bubble-bearing magmas［J］. J. Volcanol. Geotherm. Res.，1998，87：15-28.

［8］ SUNG W J，SEUNG H L，JAE R Y. Rheological behavior of polydispersed bubble suspensions in shear flows［J］. Polymer Engineering and Science，2010，50（1）：128-137

［9］ PAL R. rheological constitutive equation for bubbly suspensions［J］. Ind. Eng. Chem. Res.，2004，43：5372-5379.

［10］ FRANKEL N A，ACRIVOS A. The constitutive equation for a dilute emulsion［J］. J. Fluid Mech.，1970，44：65-78.

［11］ YUN M L，DONG J S，JAE R Y. Rheological behavior of dilute bubble suspensions in polyol［J］. Korea-Australia Rheology Journal，2004，16：47-54.

［12］ HIRT C W，NICHOLS B D. Volume of fluid（VOF）method for the dynamics of free boundaries［J］. Journal of Computational Physics， 1981，39：201-205.

［13］ RIDER W J，KOTHE D B. Reconstructing volume tracking［J］. Journal of Computational Physics，1997，141：112-152.

［14］ 馬斗. 單孔氣泡動力學行為的 VOF 數(shù)值模擬［D］. 天津：天津大學化工學院，2009.

［15］ 黃萌. 電場作用下油水乳液中水滴的運動和形變特性研究［D］. 西安：西安交通大學，2015.

［16］ 李宜樂，劉應中. 繆國平. 二階精度的VOF自由面追蹤方法及其應用［J］. 船舶力學，1999，3（1）：44-52.

Numerical study on influence of dispersed bubbles on liquid-phase apparent viscosity in two-dimensional parallel plate

XU Lei，PANG Mingjun
（School of Mechanical Engineering，Changzhou University，Changzhou 213164，Jiangsu，China）

Bubble liquid exists widely in various industrial fields，so the detailed understanding of physical properties of bubbly liquid is of great significance for improving the quality of products and for strengthening industrial processes. In this paper，in order to understand the influence of bubbles on the liquid-phase apparent viscosity and to explain the changing reason of the liquid-phase apparent viscosity，based on two-dimensional parallel plate model，the effect of large bubbles on the liquid-phase apparent viscosity was deeply investigated with the volume of fluid （VOF） method combined with the dynamic mesh. The effect of bubbles on the liquid-phase apparent viscosity was numerically studied in detail under the conditions with the same void fraction but different capillary numbers. The present study showed that the capillary number and the bubble volume fraction has the great effect on liquid-phase apparent viscosity. For the same volume fraction，when the capillary number is relatively large （Ca＞1），the bubble injection causes the decrease of the relative viscosity of the liquid phase，and the larger the volume fraction is. The relative viscosity of the liquid phase decreases more sharper. In addition，under the same volume fraction，when the capillary number is relatively small （Ca＜1），the relative viscosity of the liquid phase increases due to the addition of bubbles，and the higher the volume fraction is. The relative viscosity of the liquid phase increases more obviously.

bubbly liquid； dispersed bubble；relative viscosity；volume of fluid（VOF）；moving grid； two-dimensional parallel plate

O 37

A

1000-6613（2016）10-3065-07

10.16085/j.issn.1000-6613.2016.10.008

2016-03-03；修改稿日期：2016-03-31。

國家自然科學基金（51376026）及江蘇省青藍工程項目。

徐磊 （1990—），碩士研究生，研究方向為泡狀流流變特性。聯(lián)系人：龐明軍，博士，副教授，研究方向為氣液兩相流動與湍流減阻。E-mail pangmj@cczu.edu.cn。