姜忠濤，韓蕊，李帥

(哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001)

氣泡在工程領(lǐng)域得到了人們廣泛的關(guān)注，用以探測石油的高壓氣槍、軍事中造成艦船毀傷的水下爆炸氣泡、螺旋槳產(chǎn)生的空泡，這一切都與氣泡的動(dòng)力學(xué)特性有著密不可分的關(guān)系。從 Rayleigh(1917)［1］對球狀脈動(dòng)氣泡的理論研究開始，人們開展了大量有關(guān)氣泡動(dòng)力學(xué)的研究工作，目前國內(nèi)外學(xué)者在數(shù)值和實(shí)驗(yàn)研究方面也取得了豐碩的成果［2-6］，而在數(shù)值研究方面，基于勢流理論的邊界元方法由于其高效、高精度的優(yōu)勢得到了廣泛應(yīng)用。對于單氣泡的研究，學(xué)者模擬了氣泡與自由液面［2，7-9］、剛性壁面［10］等不同邊界的相互作用，并且對射流后的環(huán)狀氣泡也進(jìn)行了數(shù)值模擬［11-12］;而隨著研究的深入，很多學(xué)者對2個(gè)氣泡與自由液面［13-14］、剛性壁面［15］的工況也進(jìn)行了研究，F(xiàn)ong 等人(2008)［16］通過電火花實(shí)驗(yàn)對氣泡相位差、氣泡間距聯(lián)合作用下2個(gè)氣泡相互作用的不同現(xiàn)象進(jìn)行了圖表總結(jié)，同時(shí)采用邊界元方法只對3個(gè)氣泡相互作用的幾個(gè)典型工況進(jìn)行了數(shù)值模擬。研究3個(gè)氣泡相互作用的公開文獻(xiàn)還較少，這主要是由于氣泡運(yùn)動(dòng)與氣泡間距、初始條件、排列方式、相位等多因素有關(guān)［16］，其中較基礎(chǔ)的研究是垂向布置的3個(gè)同相氣泡間的相互作用。

本文基于勢流理論，建立了邊界元數(shù)值模型，對3個(gè)氣泡的相互作用進(jìn)行數(shù)值模擬。首先，將數(shù)值結(jié)果與球狀脈動(dòng)氣泡解析解進(jìn)行對比，驗(yàn)證了數(shù)值模型的有效性和收斂性，然后研究垂向布置的3個(gè)同相氣泡的耦合特性，主要分析氣泡間距對氣泡運(yùn)動(dòng)的影響，分別探討了等間距和不等間距布置氣泡工況中氣泡的運(yùn)動(dòng)規(guī)律和機(jī)理，旨為研究多氣泡相互作用的研究提供一定的參考。

1 邊界元基本理論

本文基于勢流理論，建立軸對稱數(shù)值模型模擬垂向布置的3個(gè)同相氣泡間的相互作用，由于氣泡運(yùn)動(dòng)過程中，氣泡表面速度維持在比水中聲速小一些的量級上，且雷諾數(shù)高、持續(xù)時(shí)間短，故可忽略可壓縮性和粘性耗散。因此，認(rèn)為流體不可壓縮、無粘、無旋，速度勢在流域內(nèi)滿足拉普拉斯方程以及邊界積分方程:

式中:p是邊界上控制點(diǎn);ω-為在該點(diǎn)觀察流場的立體角，此處=2π;q則為是邊界上積分點(diǎn);G(p，為自由空間格林函數(shù);n是指向流場外(氣泡內(nèi)部)的法向量，法向?qū)?shù)定義為為流體域的所有邊界，在本文中所有邊界由3個(gè)氣泡的表面組成。

由于氣體和水的密度巨大差異，所以在本文中假設(shè)氣泡內(nèi)氣體運(yùn)動(dòng)對壓力的影響可忽略不計(jì);同時(shí)由于氣泡運(yùn)動(dòng)的瞬態(tài)特性，故可忽略熱交換效應(yīng)，基于絕熱假設(shè)，氣體壓力僅和氣泡初始狀態(tài)及其體積有關(guān)，則每個(gè)氣泡外表面水的壓力均滿足［12］:

式中:Pc為氣泡內(nèi)可冷凝氣體的飽和蒸汽壓;下標(biāo)“0”表示氣泡初始狀態(tài)(初始壓力和體積);?為氣體的比熱，和氣體組成有關(guān)，本文中 ? 取 1.25［6，17?。

氣泡表面的動(dòng)力學(xué)邊界條件為

式中:P∞為靜水壓力，ρ為流體密度，g為重力加速度，z代表位置矢量的垂直分量。

研究過程中，采用氣泡最大半徑Rm、ΔP=P∞-Pc分別作為長度、壓力的特征量，相應(yīng)地以Rm(ρ/ΔP)1/2、Rm(ΔP/ρ)1/2以及 (ρgRm/ΔP)1/2對時(shí)間、速度勢和浮力進(jìn)行無量綱化，則得到式(4)的無量綱形式為［18］

式中:ε=P0/ΔP為氣泡強(qiáng)度參數(shù)為浮力參數(shù)，本文中κ=0。同時(shí)，對于本文中研究的3個(gè)氣泡相互作用，如圖1垂向布置3個(gè)氣泡，則分別定義d1=D1/Rm和d2=D2/Rm為氣泡1、2和氣泡2、3之間的間距參數(shù)，D1和D2為相應(yīng)的氣泡中心間距，在等間距工況中間距參數(shù)滿足d1=d2=d。另外，下文中出現(xiàn)的其他無量綱變量還包括時(shí)間T、氣泡半徑R、射流速度vjet和氣泡體積V。圖1中參數(shù)R0為氣泡的初始半徑，可由下式確定［19］:

氣泡的運(yùn)動(dòng)學(xué)邊界條件為［11］

數(shù)值求解方程(2)的具體細(xì)節(jié)可參考文獻(xiàn)［19］，根據(jù)式(2)求得的速度，聯(lián)立式(5)和式(7)，則可對氣泡表面的速度勢及位置進(jìn)行更新。計(jì)算過程中采用二階龍格-庫塔沿時(shí)域向前推進(jìn)，且為維持計(jì)算過程的穩(wěn)定，每一時(shí)間步長均需嚴(yán)格控制，使得每個(gè)節(jié)點(diǎn)的速度勢在每一步的改變量不會(huì)超過Δφ(Δφ為某一常數(shù))，其取值將在下文進(jìn)行討論說明。

圖1 氣泡布置及間距定義Fig.1 Bubble arrangement and the definition of inter-bubble distance

2 數(shù)值結(jié)果

2.1 數(shù)值模型有效性驗(yàn)證

本文建立邊界元模型研究多氣泡的耦合特性，在此之前需要驗(yàn)證數(shù)值模型的有效性。在本文中，首先采用Rayleigh-Plesset方程(R-P方程)描述無重力自由場中的球形氣泡的運(yùn)動(dòng)，R-P方程的無量綱形式為［1］

選取工況為 ε=60，R0=0.178 8，κ=0，采用四階龍格-庫塔求解，可得到球形氣泡半徑時(shí)歷變化曲線。同時(shí)，用本文的邊界元模型對相同工況下球形氣泡的運(yùn)動(dòng)進(jìn)行計(jì)算，將得到的結(jié)果與R-P方程對比，如圖2所示。圖2(a)和圖2(b)分別給出了節(jié)點(diǎn)數(shù)和ΔΦ的改變對計(jì)算結(jié)果的影響，觀察發(fā)現(xiàn)數(shù)值結(jié)果與R-P方程解析解均吻合良好。在不考慮重力作用的自由場中，氣泡的最大半徑理論值為1，分別列出不同節(jié)點(diǎn)數(shù)和Δφ的計(jì)算結(jié)果、相對誤差及計(jì)算時(shí)長，見表1、表2。

由表1可知，當(dāng)氣泡節(jié)點(diǎn)數(shù)增加至30以上時(shí)，最大半徑Rm的誤差已減小至0.1%以下，隨著節(jié)點(diǎn)數(shù)的增加，計(jì)算時(shí)長增加;由表2可知，當(dāng)Δφ≤0.03(節(jié)點(diǎn)數(shù)為50)時(shí)，計(jì)算誤差都能保持0.05%以內(nèi)，隨著Δφ的減小，計(jì)算精度提高并不明顯，但計(jì)算時(shí)長明顯增加。

圖2 氣泡半徑時(shí)歷變化曲線數(shù)值結(jié)果與解析解對比Fig.2 Comparison of the bubble radius history of the BEM model with the analytical solution

表1 不同節(jié)點(diǎn)數(shù)的計(jì)算精度及效率Table 1 Effect of nodes number on the calculation accuracy and efficiency

表2 不同Δφ的計(jì)算精度及效率Table 2 Effect of Δφ on the calculation accuracy and efficiency

通過圖2曲線對比和表1、表2的數(shù)據(jù)可驗(yàn)證本文數(shù)值模型的有效性，同時(shí)，數(shù)值模型的收斂性也可以得到保證。在下文數(shù)值算例中，選擇節(jié)點(diǎn)數(shù)50和Δφ=0.02以平衡計(jì)算精度和計(jì)算效率(計(jì)算時(shí)長)。

2.2 數(shù)值算例

數(shù)值模擬3個(gè)同相氣泡間的耦合特性，3個(gè)初始條件相同的氣泡垂向排列，初始參數(shù)均為ε=60，R0=0.1788，κ =0，每個(gè)氣泡表面采用 50 節(jié)點(diǎn)，取Δφ=0.02，分別對等間距和不等間距的情況進(jìn)行計(jì)算分析。

圖3 不同d下，射流完成時(shí)氣泡形態(tài)Fig.3 Bubble profiles upon the impact with different d

首先研究等間距情況，即d1=d2=d。由于氣泡間距過小時(shí)，氣泡會(huì)發(fā)生融合的現(xiàn)象，而間距過大時(shí)則氣泡間影響不明顯，故本文中氣泡間距的變化范圍取為［1.8，5.0］。當(dāng)氣泡間距d發(fā)生變化時(shí)，氣泡射流完成時(shí)的形態(tài)也會(huì)有所不同，如圖3所示。觀察發(fā)現(xiàn)，當(dāng)d=1.8時(shí)，3個(gè)氣泡的相鄰位置在很大范圍內(nèi)均發(fā)生扁平的現(xiàn)象，氣泡2更呈現(xiàn)出“圓柱狀”，且其中間部分出現(xiàn)向內(nèi)的凹陷;d=2.0工況中的氣泡形態(tài)與上一工況相似，只是扁平范圍較小且氣泡2較“瘦長”;當(dāng)氣泡間距d繼續(xù)增大時(shí)，可以發(fā)現(xiàn)氣泡1和3在射流完成時(shí)的體積減小，且被射流沖擊區(qū)域已不再扁平，而氣泡2也已呈現(xiàn)“橢球狀”并逐漸趨于球形，同時(shí)中間位置也不再出現(xiàn)凹陷。為了更準(zhǔn)確地分析氣泡形態(tài)的變化，計(jì)算不同d對應(yīng)的氣泡1射流速度時(shí)歷曲線、氣泡1體積變化和氣泡2體積變化，如圖4～6所示。

由于氣泡1、3的運(yùn)動(dòng)相同，其形態(tài)以氣泡2為中心對稱，所以只給出氣泡1的計(jì)算結(jié)果。由圖4的射流速度時(shí)歷曲線可了解到氣泡的運(yùn)動(dòng)過程，圖中實(shí)線為射流點(diǎn)(氣泡1頂點(diǎn))速度、虛線為被射流沖擊點(diǎn)(氣泡1最低點(diǎn)，與氣泡2相鄰點(diǎn))速度，氣泡首先快速膨脹，隨著內(nèi)部壓力的減小，膨脹速度變慢、膨脹運(yùn)動(dòng)減緩，可以觀察到膨脹階段氣泡基本保持球形;當(dāng)膨脹到某個(gè)時(shí)刻，膨脹停止，緊接著氣泡快速收縮，由于周圍另一氣泡的影響，氣泡上、下兩部分收縮程度不同，隨著收縮的進(jìn)行，氣泡1形成向下的射流，同時(shí)可觀察到，氣泡收縮速度在其后期變慢，這主要是由于氣泡收縮導(dǎo)致內(nèi)壓增大，從而收縮運(yùn)動(dòng)受到阻礙。觀察圖4射流速度受d影響的變化規(guī)律:d不同，射流點(diǎn)速度在坍塌階段后期有所不同，隨著d的增大，氣泡射流點(diǎn)速度在后期增加越快、射流完成時(shí)速度越大;隨著d的增大，氣泡被射流沖擊點(diǎn)的膨脹速度減小得慢，坍塌階段收縮速度快(如圖5所示，氣泡1達(dá)到的最大體積變大，但體積減小速度快，從而射流完成時(shí)間早)，但受內(nèi)壓影響較大，導(dǎo)致坍塌后期收縮速度減小，尤其d∈［3.0，5.0］范圍內(nèi)收縮速度減小十分明顯。由以上敘述可知，當(dāng)氣泡周圍存在另一氣泡時(shí)，其運(yùn)動(dòng)將受到抑制，從而氣泡能到達(dá)的最大體積變小，但其收縮速度受到的抑制導(dǎo)致氣泡射流完成得晚。以同樣的原因分析氣泡2的體積變化，觀察圖6發(fā)現(xiàn)，d較小時(shí)，氣泡2受到來自上、下方氣泡的抑制，其膨脹速度慢，但氣泡2在此作用下形狀會(huì)發(fā)生變化，呈現(xiàn)出被拉長的“橢球狀”或上下端扁平的“圓柱狀”，故膨脹階段氣泡2能夠達(dá)到的最大體積增大，但d對接下來的收縮速度影響不大。由此可知，多氣泡間存在抑制作用，從而影響其膨脹和收縮速度，這種速度的不同將導(dǎo)致氣泡形態(tài)發(fā)生變化。

接下來研究不等間距情況，即d1≠d2，且由于本文設(shè)置初始參數(shù)下氣泡1、3形態(tài)、運(yùn)動(dòng)的對稱性，故只需研究d1≤d2的情況。為避免融合現(xiàn)象的發(fā)生，不等間距情況的研究中d1∈ [1.8，5]且d1+d2=10，在此間距組合下，由以上分析可知?dú)馀?的運(yùn)動(dòng)受到抑制，計(jì)算在氣泡3完成射流時(shí)刻停止，故圖7給出氣泡3的射流速度時(shí)歷曲線，同時(shí)給出氣泡1和氣泡3體積變化以及氣泡2體積變化，如圖8、9所示。圖7、圖8中的變化規(guī)律與等間距分析的結(jié)果相同，在此不再贅述:仔細(xì)觀察圖7可發(fā)現(xiàn)被射流沖擊點(diǎn)(氣泡3頂點(diǎn))的速度在坍塌階段后期減小十分明顯，與等間距情況中d∈［3.0，5.0］的變化類似，說明氣泡3已經(jīng)開始發(fā)生回彈(圖8中氣泡3體積變化曲線已經(jīng)開始上升);觀察圖8可知，隨著d1，d2差距的減小，氣泡1和氣泡3體積變化曲線均向d1=d2=5.0的氣泡1體積變化曲線移動(dòng)，說明了該數(shù)值模型和變化規(guī)律的收斂性。圖9給出了氣泡2的體積變化曲線，可以發(fā)現(xiàn)d1，d2的變化對氣泡2達(dá)到的最大體積影響不大，而隨著d1的減小和d2的增大，氣泡2的膨脹、收縮速度都增大。

圖4 氣泡1射流速度時(shí)歷曲線Fig.4 Time-history of the jet velocity of bubble 1

圖5 氣泡1體積變化Fig.5 Time-history of the volume of bubble 1

圖6 氣泡2體積變化Fig.6 Time-history of volume of bubble 2

圖7 氣泡3射流速度時(shí)歷曲線Fig.7 Time-history of the jet velocity of bubble 3

圖8 氣泡1和氣泡3體積變化Fig.8 Time-history of the volume of bubble 1 and bubble 3

圖9 氣泡2體積Fig.9 Time-history of the volume of bubble 2

3 結(jié)論

1)多氣泡間存在抑制作用，垂向等間距布置的3個(gè)氣泡，上、下兩端氣泡會(huì)出現(xiàn)射流現(xiàn)象，中間氣泡則在兩端氣泡的影響下呈現(xiàn)出“圓柱狀”或“橢球狀”，若間距足夠大則趨于球形脈動(dòng)。

2)3個(gè)同相氣泡垂向等間距布置，隨著d的增大，上、下兩端氣泡能達(dá)到的最大體積越大，但完成射流越早，而隨著d的增大中間氣泡能夠達(dá)到的最大體積越小。

3)在不等間距布置氣泡的情況下，最下方氣泡(與中間氣泡間隔遠(yuǎn))首先發(fā)生射流，隨著中間氣泡不斷向下移動(dòng)，下方氣泡的膨脹、收縮運(yùn)動(dòng)均受到阻礙，而中間氣泡的膨脹、收縮速度都變大，但其最大體積變化不明顯。

［1］RAYLEIGH L.VIII.On the pressure developed in a liquid during the collapse of a spherical cavity［J］.Philosophical Magazine Series 6，1917，34(200):94-98.

［2］BLAKE J R，GIBSON D C.Growth and collapse of a vapour cavity near a free surface［J］.Journal of Fluid Mechanics，1981，111:123-140.

［3］WANG Q X，YEO K S，KHOO B C，et al.Vortex ring modelling of toroidal bubbles［J］.Theoretical and Computational Fluid Dynamics，2005，19(5):303-317.

［4］KLASEBOER E，HUNG K C，WANG C，et al.Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/rigid structure［J］.Journal of Fluid Mechanics，2005，537:387-413.

［5］ZHANG Aman，YANG Wenshan，HUANG Chao，et al.Numerical simulation of column charge underwater explosion based on SPH and BEM combination［J］.Computers＆ Fluids，2013，71:169-178.

［6］ZHANG Aman，WANG Shiping，HUANG Chao，et al.Influences of initial and boundary conditions on underwater explosion bubble dynamics［J］.European Journal of Mechanics-B/Fluids，2013，42:69-91.

［7］BLAKE J R，TAIB B B，DOHERTY G.Transient cavities near boundaries.Part 2.Free surface［J］.Journal of Fluid Mechanics，1987，181:197-212.

［8］WANG Q X，YEO K S，KHOO B C，et al.Strong interaction between a buoyancy bubble and a free surface［J］.Theoretical and Computational Fluid Dynamics，1996，8(1):73-88.

［9］李帥，張阿漫，王詩平.氣泡引起的皇冠型水冢實(shí)驗(yàn)與數(shù)值研究［J］.物理學(xué)報(bào)，2013，62(19):194703.LI Shuai，ZHANG Aman，WANG Shiping Experimental and numerical studies on“crown”spike generated by a bubble near free-surface［J］.Acta Physica Sinica，2013，62(19):194703.

［10］LI Shuai，LI Yunbo，ZHANG Aman.Numerical analysis of the bubble jet impact on a rigid wall［J］.Applied Ocean Research，2015，50:227-236.

［11］BEST J P.The formation of toroidal bubbles upon the collapse of transient cavities［J］.Journal of Fluid Mechanics，1993，251:79-107.

［12］WANG Q X，YEO K S，KHOO B C，et al.Nonlinear interaction between gas bubble and free surface［J］.Computers＆ Fluids，1996，25(7):607-628.

［13］ROBINSON P B，BLAKE J R，KODAMA T，et al.Interaction of cavitation bubbles with a free surface［J］.Journal of Applied Physics，2001，89(12):8225-8237.

［14］HAN Rui，ZHANG Aman，LI Shuai.Three-dimensional numerical simulation of crown spike due to coupling effect between bubbles and free surface［J］.Chinese Physics B，2014，23(3):034703.

［15］BLAKE J R，ROBINSON P B，SHIMA A，et al.Interaction of two cavitation bubbles with a rigid boundary［J］.Journal of Fluid Mechanics，1993，255:707-721.

［16］FONG S W，ADHIKARI D，KLASEBOER E，et al.Interactions of multiple spark-generated bubbles with phase differences［J］.Experiments in Fluids，2009，46(4):705-724.

［17］李章銳，宗智，董婧，等.兩個(gè)氣泡相互作用的某些動(dòng)力特性研究［J］.船舶力學(xué)，2012，16(7):717-729.LI Zhangrui，ZONG Zhi，DONG Jing，et al.Some dynamic characteristics of the interactions of two bubbles［J］.Journal of Ship Mechanics，2012，16(7):717-729.

［18］WANG C，KHOO B C.An indirect boundary element method for three-dimensional explosion bubbles［J］.Journal of Computational Physics，2004，194(2):451-480.

［19］KLASEBOER E，HUNG K C，WANG C，et al.Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/rigid structure［J］.Journal of Fluid Mechanics，2005，537:387-413.