多氣泡耦合作用下非球狀特性數(shù)值研究

姜忠濤，韓蕊，李帥

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

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

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

1 邊界元基本理論

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

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

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

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

氣泡表面的動力學邊界條件為

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

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

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

氣泡的運動學邊界條件為［11］

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

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

2 數(shù)值結果

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

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

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

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

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

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

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

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

2.2 數(shù)值算例

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

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

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

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

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

圖4 氣泡1射流速度時歷曲線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射流速度時歷曲線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 結論

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

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

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

［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ù)值研究［J］.物理學報，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］李章銳，宗智，董婧，等.兩個氣泡相互作用的某些動力特性研究［J］.船舶力學，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.