梁坤峰，王全海，阮春蕾，王 林

(河南科技大學(xué)a.車輛與動力工程學(xué)院;b.數(shù)學(xué)與統(tǒng)計學(xué)院;c.規(guī)劃與建筑工程學(xué)院，河南洛陽471003)

0 引言

近年來，人們在結(jié)晶、流體力學(xué)和外界因素干預(yù)之間相互關(guān)聯(lián)方面的認(rèn)識取得了重大的進步，已開始積極尋找各種外界因素動態(tài)干預(yù)結(jié)晶過程與環(huán)境［1－5］。外加超聲強化水或溶液結(jié)晶的研究越來越受到科研工作者的重視，雖然超聲強化水或水溶液結(jié)晶的原因仍不十分清楚，但普遍認(rèn)為超聲作用于兩相或多相體系產(chǎn)生的空化效應(yīng)，即液體中存在或生成的氣體微泡和空穴微泡，在超聲場作用下振動、縮放和崩潰等一系列動力學(xué)過程促進了晶體的形成［4］。由于微氣泡廣泛存在于化工、動力設(shè)備、熱能、醫(yī)藥等領(lǐng)域，對微氣泡運動規(guī)律的研究對實際生產(chǎn)、流體力學(xué)及多相流動的發(fā)展具有重要意義。

目前，對氣泡的上浮速度、氣泡分布及超聲應(yīng)用等問題的研究已受到國內(nèi)外學(xué)者廣泛的關(guān)注。文獻［6］測定了單個氣泡在水及不同黏度甘油水溶液中的形變特性及上升終速度;文獻［7］研究了氣泡半徑變化和上升運動的耦合影響，建立了運動方程，并用數(shù)值計算得出了氣泡半徑變化對上升運動的影響;文獻［8］從力學(xué)平衡方程的角度，對氣泡上浮過程中所受的各種主要力學(xué)影響因子進行了考慮，并詳細(xì)討論了黏滯系數(shù)與氣泡運動阻力系數(shù)對氣泡運動的重要影響;文獻［9］研究了單個氣泡在遠(yuǎn)程尾流中的運動特性及影響要素，重點考慮了上浮速度與傳質(zhì)速率相互影響的耦合因素。然而，上述研究在考慮氣泡上浮過程中所受的主要力學(xué)影響因子，氣泡半徑變化對氣泡上升運動特性影響以及外加超聲場作用下微氣泡運動特性等方面的研究尚不完善。

本文從力學(xué)平衡角度出發(fā)，建立了水中單個氣泡的運動平衡方程，研究了各種力學(xué)影響因子對單個氣泡上浮過程運動特性的影響;分析了不同粒徑氣泡時上浮速度的變化規(guī)律;進而重點探討了氣泡在外加超聲波作用下，氣泡粒徑的變化規(guī)律以及氣泡的運動特性。

1 模擬理論與方法

氣泡運動過程由其受力情況決定，為了研究方便，做出以下假設(shè):

(Ⅰ)氣泡在運動過程中保持球形狀態(tài)。

(Ⅱ)氣泡內(nèi)氣體溫度不變。

(Ⅲ)液體處于靜止?fàn)顟B(tài)，無流動。

根據(jù)假設(shè)，分析氣泡的各個受力，建立單個氣泡運動平衡方程，如下:

其中，F(xiàn)g為氣泡的重力;Fb為氣泡受到的浮力;Fd為運動黏性阻力;FA為虛擬質(zhì)量力;FB為Basset 力［10］。

將各力的表達(dá)式代入平衡方程式(1)，經(jīng)過整理化簡可得:

由氣泡運動微分方程式(2)可知:式中關(guān)于Basset 力的積分式包含一奇異端點，對其求解需要妥善處理該廣義積分項。文獻［11］基于廣義積分的極限審斂法，判定該廣義積分項必定收斂，提出小區(qū)間［t－△t］內(nèi)近似公式，消除奇異點，直接由復(fù)化梯形公式求解積分項，如式(3)所示。

為提高數(shù)值計算結(jié)果的精度，采用先預(yù)測再校正的方法即改進的Eular 法［12－13］，通過自適應(yīng)變步長聯(lián)立求解方程。即當(dāng)解的變化較慢時采用較大的計算步長，從而使計算速度加快;當(dāng)方程的解變化較快時，積分步長適當(dāng)變小，從而使得計算的精度提高。模型求解中，一些關(guān)鍵參數(shù)的取值如表1 所示。

表1 各參數(shù)取值

2 計算結(jié)果及分析

2.1 力學(xué)因子影響

圖1 為氣泡在力的影響下上浮運動速度隨時間的變化曲線，由圖1 可以看出:隨著時間的增加速度逐漸增加，運動到某一時刻速度趨于穩(wěn)定。在圖1 中，曲線a 表示考慮虛擬質(zhì)量力和Basset 力時的速度曲線，曲線b 表示不考慮Basset 力時的速度曲線，對比這兩條曲線，發(fā)現(xiàn)不考慮Basset 力時，氣泡運動速度比考慮Basset 力時速度大很多，表明Basset 力不僅持續(xù)影響氣泡的運動，而且起主要作用。而在不考慮虛擬質(zhì)量力和Basset 力及不考慮虛擬質(zhì)量力時，氣泡上浮過程的初始階段，氣泡能從零突變到穩(wěn)態(tài)速度，不考慮Basset 力與不考慮虛擬質(zhì)量力及Basset 力時的曲線在3 ms 時刻后基本重合，表明虛擬質(zhì)量力及Basset 力對氣泡的運動有影響。此外，為清晰展示初始階段不考慮虛擬質(zhì)量力和Basset 力及不考慮虛擬質(zhì)量力時氣泡上浮速度突變的規(guī)律，圖2 給出了氣泡上浮初始階段各個力的影響曲線放大視圖;由圖2 可知:氣泡運動初期，速度變化非常劇烈，具有顯著的震蕩特征，但很快趨于穩(wěn)定，時間為0.01 ms，由此進一步表明虛擬質(zhì)量力和Basset 力需要同時考慮。

圖1 氣泡上浮過程各個力的影響曲線

圖2 氣泡上浮初始階段各個力的影響曲線

2.2 氣泡粒徑影響

圖3 為粒徑分別取30 μm、50 μm、100 μm、200 μm、300 μm 時，氣泡上浮過程的速度隨時間的變化曲線。由圖3 可見:不同粒徑，隨著時間的增加氣泡速度均呈遞增趨勢，且初始時刻速度增長較快，在氣泡上浮過程后期其速度大小趨于穩(wěn)定;同一時刻，氣泡的粒徑越大其上浮速度越大，反之，氣泡的粒徑越小其上浮速度越小;此外，當(dāng)氣泡粒徑小于50 μm 時，隨著時間的增加，氣泡上浮速度增加趨勢不明顯。圖3 曲線表明:氣泡的上浮速度與氣泡自身的粒徑有直接關(guān)系，因為粒徑大的氣泡受到液體的浮力大，所以氣泡運動速度就快，而對于粒徑小的氣泡，即氣泡粒徑小于50 μm 時，如要獲得較大的速度，則須借助外力作用。

圖3 定粒徑氣泡速度變化圖

表2 為氣泡粒徑分別取30 μm、50 μm、100 μm、200 μm、300 μm 時，氣泡運動末速度的實驗值與模擬值的對比表。模擬值可以根據(jù)方程(3)計算得出，實驗值則依據(jù)文獻［4］可得，由表2 可以看出:實驗值與模擬值的變化規(guī)律一致，即氣泡的粒徑越大其上浮速度越大，且氣泡上浮過程中，模擬值與實驗值誤差均小于8.0%，表明本文所提出的求解方法可行。

2.3 超聲波影響

當(dāng)外加超聲Pasin［2πf(t－x/c)］作用于氣泡時，氣泡受到激勵聲波的拉伸和壓縮作用產(chǎn)生振動。假設(shè)液體是理想的、不可壓縮的連續(xù)介質(zhì)，黏滯系數(shù)為常數(shù)，并設(shè)激勵聲場頻率和強度恒定，聲波波長遠(yuǎn)大于空化氣泡半徑，氣泡在運動中內(nèi)部蒸汽壓保持常數(shù)，與周圍液體無熱量交換，可由流體Navier-Stokes 方程、狀態(tài)方程、無限處和氣泡壁處的邊界條件而建立RP 方程［14］。

表2 實驗與模擬結(jié)果對比

式中，r0為氣泡初始半徑;r 為氣泡運動半徑;ρ 為氣泡外液體的密度;σ 為氣泡外液體的表面張力;p0為氣泡外大氣壓;pv為水蒸氣壓;pa為超聲激勵聲壓;μ 為液體的黏度系數(shù);γ 為等壓比熱容與等容比熱容的比值。各參數(shù)取值如表3 所示。

表3 參數(shù)取值

在超聲波作用下，氣泡既受重力、浮力、黏性阻力、虛擬質(zhì)量力和Basset 力的影響，還會受到聲壓梯度力的作用，位于x 處聲壓梯度力可由式(5)計算［15］:

此時，氣泡的運動平衡方程變?yōu)?

圖4 為聲波作用下氣泡速度變化曲線。設(shè)初始時刻氣泡處于超聲波1/8 波長處，即x =1/8λ，采用Runge-Kutta 法數(shù)值求解方程(4)，得到氣泡粒徑隨時間t 的變化關(guān)系，如圖4a 所示，表明在超聲波作用下氣泡的粒徑呈震蕩波動變化。根據(jù)新的平衡方程式(6)，采用改進的Eular 法模擬求解，得到超聲波作用下氣泡上浮速度曲線，其結(jié)果如圖4b 所示，其中，實線為考慮聲壓梯度力時氣泡速度變化曲線，其縱坐標(biāo)對應(yīng)左側(cè)的縱坐標(biāo);虛線為不考慮聲壓梯度力時氣泡的速度變化曲線，其縱坐標(biāo)對應(yīng)右側(cè)的縱坐標(biāo)。由圖4b 可知:超聲波作用下，不考慮聲壓梯度力時，氣泡速度變化整體上呈遞增趨勢，但其速度值量級較小;當(dāng)考慮聲壓梯度力時，氣泡速度呈近似正弦波動狀態(tài)，且其振動的振幅逐漸增大，速度數(shù)量級比不考慮聲壓梯度力時的大很多;表明超聲波作用下，氣泡受到聲壓梯度力作用，并對氣泡的運動產(chǎn)生很大的影響，研究氣泡運動時必須考慮超聲引起的聲壓梯度力。

圖4 聲波作用下氣泡速度變化曲線

圖5 為初始位置x=1/4λ、1/2λ、3/4λ、λ 時，氣泡速度隨時間變化的曲線圖。由圖5 可知:不同的初始位置x 處，其速度變化規(guī)律不同，其原因為在外加超聲作用時，氣泡受聲壓梯度力FP的大小、方向與氣泡初始時刻的位置x 有關(guān)。對比x=1/4λ 與x=3/4λ 時的速度曲線可以看出:兩條曲線是對稱的，同樣，x=1/2λ 與x=λ 時的速度曲線也是對稱的，因為在x=1/4λ 與x=3/4λ 和x=1/2λ 與x=λ 處氣泡受聲壓梯度力FP的方向相反;對比x =1/4λ 與x =1/2λ 時的速度曲線可見:x =1/4λ 時的速度峰值比x=1/2λ 時速度峰值略高，因為x=1/4λ 時的聲壓梯度力FP較大。

圖6 為初始時刻x=1/4λ、1/2λ、3/4λ、λ 時，氣泡上浮運動的位移隨時間的變化曲線圖。由圖6 可見:氣泡速度的波動變化引起氣泡位移的波動;當(dāng)x =1/4λ、1/2λ 時，運動后期氣泡的位移向下，x =3/4λ、λ 時，運動后期氣泡的位移向上;x=1/4λ 與x =3/4λ時，氣泡的位移曲線對稱，同樣，在x =1/2λ與x=λ 時，氣泡的位移曲線亦呈對稱特征，均因其所受的聲壓梯度力FP的方向相反。

圖5 氣泡初始位置對速度的影響

圖6 氣泡初始位置對位移的影響

3 結(jié)論

(1)靜水中氣泡的上浮速度隨時間增加而逐漸增大，計算值與實驗值具有很好的一致性。上浮初始階段必須考慮由于加速度所引起的虛擬質(zhì)量力和Basset 力對上浮速度的影響，且在不考慮虛擬質(zhì)量力時，上浮速度具有震蕩特征;氣泡上浮后期，上浮速度趨于穩(wěn)定;在所計算的氣泡粒徑范圍(30 ～300 μm)內(nèi)，氣泡上浮速度從靜止變化到穩(wěn)定所經(jīng)歷的時間均在10 ms 量級。

(2)靜水中氣泡初始半徑越大，氣泡上浮速度越大，氣泡初始粒徑越小，氣泡上浮速度越小，但達(dá)到穩(wěn)定的時間越短;當(dāng)氣泡初始粒徑小于50 μm 時，依靠氣泡上浮的各種力學(xué)影響因子難以達(dá)到較大的上浮速度。

(3)聲壓梯度力對氣泡運動的影響強于其他力學(xué)影響因子，速度相差兩個數(shù)量級;氣泡速度的變化與超聲波發(fā)生規(guī)律一致，近似呈現(xiàn)正弦波動的變化規(guī)律，速度波動的振幅逐漸增大;氣泡的初始位置處于超聲波的不同波長處，氣泡運動的方向相異，使得氣泡具有上下兩個方向的位移變化。

［1］ Yu D Y，Liu B L，Wang B C.The Effect of Ultrasonic Waves on the Nucleation of Pure Water and Degassed Water［J］.Ultrasonics Sonochemistry，2012，19(3):459－463.

［2］ Heneghan A F，Haymet A D J. Liquid-to-crystal Heterogeneous Nucleation Bubble Accelerated Nucleation of Pure Supercooled Water［J］.Chemical Physics Letters，2003，368(1/2):177－182.

［3］ Mathieu S，Roman P，Julien A. A Theoretical Model for Ice Primary Nucleation Induced by Acoustic Cavitation［J］.Ultrasonics Sonochemistry，2010，17(1):98－105.

［4］ Thomas L，Muthupandian A，Sandra K.The Fundamentals of Power Ultrasound:A Review［J］.Acoustics Australia，2011，39(2):54－63.

［5］ 李玉剛，王程遠(yuǎn)，陳衛(wèi)鋒，等.液－液霧化液滴的粒徑分布特性［J］.河南科技大學(xué)學(xué)報:自然科學(xué)版，2010，31(3):29－32.

［6］ 郭容，蔡子琦，高正明.黏性流體中單氣泡的運動特性［J］.高校化學(xué)工程學(xué)報，2009，23(6):916－921.

［7］ Chakraborty B B，Tuteja G S.Motion of an Expanding Spherical Gas Bubble in a Viscous Liquid Under Gravity［J］.Phys Fluids A，1993，5(8):1879－1882.

［8］ 石晟瑋，王江安.水中微氣泡上浮過程的力學(xué)影響因子研究［J］.海軍工程大學(xué)學(xué)報，2008，6(3):83－87.

［9］ 田恒斗，金良安，王涌，等.考慮單氣泡運動特性的艦船尾流氣泡分布研究［J］.兵工學(xué)報，2011，9(9):1126－1131.

［10］ Zhang J，F(xiàn)an L S.On the Rise Velocity of an Interactive Bubble in Liquids［M］.Chemical Engineering Journal，2003，92:169－176.

［11］ Carrica P M，Bonetto F J，Drew D A，et al.The Interaction of Background Ocean Air Bubbles with a Surface Ship［J］.Number Meth Fluids，1998，28(4):571－600.

［12］ 李慶揚，王能超，易大義.數(shù)值分析［M］.武漢:華中科技大學(xué)出版社，2002:88－90.

［13］ Chen M.An Automatic Recognition of the Best Threshold for Image［J］. Computer Applications and Software，2006，23(4):85－86.

［14］ Stride B E.The Influence of Surface Adsorption on Microbubble Dynamics［J］.Philos Trans A Math Phys Eng Sci，2008，366(1873):2013－2015.

［15］ 袁竹林，李偉力.細(xì)微顆粒在行波和駐波聲場中運動特性數(shù)值實驗［J］.東南大學(xué)學(xué)報，2005，35(1):140－144.