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

姜忠濤，韓蕊，李帥

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

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

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

1 邊界元基本理論

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

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

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

式中:Pc為氣泡內可冷凝氣體的飽和蒸汽壓;下標“0”表示氣泡初始狀態(初始壓力和體積);?為氣體的比熱，和氣體組成有關，本文中 ? 取 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為氣泡強度參數為浮力參數，本文中κ=0。同時，對于本文中研究的3個氣泡相互作用，如圖1垂向布置3個氣泡，則分別定義d1=D1/Rm和d2=D2/Rm為氣泡1、2和氣泡2、3之間的間距參數，D1和D2為相應的氣泡中心間距，在等間距工況中間距參數滿足d1=d2=d。另外，下文中出現的其他無量綱變量還包括時間T、氣泡半徑R、射流速度vjet和氣泡體積V。圖1中參數R0為氣泡的初始半徑，可由下式確定［19］:

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

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

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

2 數值結果

2.1 數值模型有效性驗證

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

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

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

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

表1 不同節點數的計算精度及效率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的數據可驗證本文數值模型的有效性，同時，數值模型的收斂性也可以得到保證。在下文數值算例中，選擇節點數50和Δφ=0.02以平衡計算精度和計算效率(計算時長)。

2.2 數值算例

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

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

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

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

接下來研究不等間距情況，即d1≠d2，且由于本文設置初始參數下氣泡1、3形態、運動的對稱性，故只需研究d1≤d2的情況。為避免融合現象的發生，不等間距情況的研究中d1∈ [1.8，5]且d1+d2=10，在此間距組合下，由以上分析可知氣泡1的運動受到抑制，計算在氣泡3完成射流時刻停止，故圖7給出氣泡3的射流速度時歷曲線，同時給出氣泡1和氣泡3體積變化以及氣泡2體積變化，如圖8、9所示。圖7、圖8中的變化規律與等間距分析的結果相同，在此不再贅述:仔細觀察圖7可發現被射流沖擊點(氣泡3頂點)的速度在坍塌階段后期減小十分明顯，與等間距情況中d∈［3.0，5.0］的變化類似，說明氣泡3已經開始發生回彈(圖8中氣泡3體積變化曲線已經開始上升);觀察圖8可知，隨著d1，d2差距的減小，氣泡1和氣泡3體積變化曲線均向d1=d2=5.0的氣泡1體積變化曲線移動，說明了該數值模型和變化規律的收斂性。圖9給出了氣泡2的體積變化曲線，可以發現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個氣泡，上、下兩端氣泡會出現射流現象，中間氣泡則在兩端氣泡的影響下呈現出“圓柱狀”或“橢球狀”，若間距足夠大則趨于球形脈動。

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

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

［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］李帥，張阿漫，王詩平.氣泡引起的皇冠型水冢實驗與數值研究［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.