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

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

0 引言

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

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

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

1 模擬理論與方法

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

(Ⅰ)氣泡在運(yùn)動(dòng)過程中保持球形狀態(tài)。

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

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

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

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

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

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

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

表1 各參數(shù)取值

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

2.1 力學(xué)因子影響

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

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

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

2.2 氣泡粒徑影響

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

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

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

2.3 超聲波影響

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

表2 實(shí)驗(yàn)與模擬結(jié)果對比

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

表3 參數(shù)取值

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

此時(shí)，氣泡的運(yùn)動(dòng)平衡方程變?yōu)?

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

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

圖5 為初始位置x=1/4λ、1/2λ、3/4λ、λ 時(shí)，氣泡速度隨時(shí)間變化的曲線圖。由圖5 可知:不同的初始位置x 處，其速度變化規(guī)律不同，其原因?yàn)樵谕饧映曌饔脮r(shí)，氣泡受聲壓梯度力FP的大小、方向與氣泡初始時(shí)刻的位置x 有關(guān)。對比x=1/4λ 與x=3/4λ 時(shí)的速度曲線可以看出:兩條曲線是對稱的，同樣，x=1/2λ 與x=λ 時(shí)的速度曲線也是對稱的，因?yàn)樵趚=1/4λ 與x=3/4λ 和x=1/2λ 與x=λ 處氣泡受聲壓梯度力FP的方向相反;對比x =1/4λ 與x =1/2λ 時(shí)的速度曲線可見:x =1/4λ 時(shí)的速度峰值比x=1/2λ 時(shí)速度峰值略高，因?yàn)閤=1/4λ 時(shí)的聲壓梯度力FP較大。

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

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

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

3 結(jié)論

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

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

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

［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é)報(bào):自然科學(xué)版，2010，31(3):29－32.

［6］ 郭容，蔡子琦，高正明.黏性流體中單氣泡的運(yùn)動(dòng)特性［J］.高?；瘜W(xué)工程學(xué)報(bào)，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é)報(bào)，2008，6(3):83－87.

［9］ 田恒斗，金良安，王涌，等.考慮單氣泡運(yùn)動(dòng)特性的艦船尾流氣泡分布研究［J］.兵工學(xué)報(bào)，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］ 李慶揚(yáng)，王能超，易大義.數(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ì)微顆粒在行波和駐波聲場中運(yùn)動(dòng)特性數(shù)值實(shí)驗(yàn)［J］.東南大學(xué)學(xué)報(bào)，2005，35(1):140－144.