基于Runge-Kutta法分析剛性平面附近空化氣泡的動力學

摘要：用鏡像空化泡等效替代剛性平面，得到超聲場中剛性平面附近空化泡的動力學方程.運用Runge-Kutta法數值計算了該方程，并和自由聲場中空化泡動力學方程的數值解進行對比分析．結果表明：剛性平面對空化氣泡脈動具有抑制作用；氣泡內的氣體絕熱系數和液體粘度系數越小，氣泡潰滅速度越大，對平面的作用也越大．研究結果有利于認識空化氣泡對剛性壁面的作用機制．

關鍵詞：Runge-Kutta法；鏡像法；剛性平面；空化氣泡

中圖分類號：O426.3 文獻標志碼：A

Analysis of Cavitation Bubbles Dynamics

Near Rigid Planes by Runge-Kutta Method

XUAN Li-wen， LIANG Jin-fu

（School of Physics and Electronic Science， Guizhou Normal University， Guiyang 550001， China）

Abstract：The dynamic equation of the cavitation bubble near the rigid plane in ultrasonic field is obtained by replacing the rigid plane with the mirror image cavitation bubble. The revised equation is numerically calculated through the Runge-Kutta method and compared with the numerical solution of the bubble’s equation in the free sound field. The results show that the rigid plane can inhibit the pulsation of cavitation bubble. The smaller the adiabatic coefficient of the gas in the bubble and the viscosity coefficient of the liquid， the greater the bubble collapse velocity and the greater the effect on the plane. These results are helpful to understand the action mechanism of cavitation bubbles on rigid walls.

Key words：Runge-Kutta method， mirror image method; rigid plane; cavitation bubble

非線性二階常微分方程在物理、化學、生物等學科中不斷出現．但由于它具有非線性，人們不能利用線性疊加原理對它進行解析求解．隨著計算機技術的發展，非線性微分方程的數值求解逐漸取代了解析求解，有效地解決了相關問題．聲空化是物理學中研究的一個熱點問題^［1］，聲空化氣泡動力學方程是一個非線性二階常微分方程^［2-3］，求精確解非常困難．本文運用四階Runge-Kutta法^［4］數值計算剛性平面附近空化氣泡的動力學方程，通過分析計算結果揭示空化氣泡對剛性壁面的作用機制．

我們知道，超聲波是一種機械波，當它通過液體時，液體內部的壓力會發生漲落變化，出現低于靜態壓力的負壓區域．當超聲波強度超過一定閾值時，液體中負壓區域會出現結構斷裂，形成缺陷．液體中的微小氣泡（空化核）在缺陷處會成長為肉眼可見的、微米量級的氣泡，這種現象稱為聲空化^［3］，形成的氣泡叫做空化泡．氣泡在聲場的作用下，在每一個聲周期內，其徑向運動表現為先緩慢膨脹，然后急劇塌縮，最后反彈的非線性特征．氣泡在坍縮期間，其最大半徑可以達到最小半徑的一百倍，氣泡的最大體積是最小體積的一百萬倍，因此，空化泡均具有異常高的聚能能力^［3］．空化泡塌縮至最小半徑前后，其內部會產生極端的高溫高壓，有時會導致光輻射，即聲致發光^［5］．空化氣泡內部的高溫高壓已被應用到超聲清洗^［6］，納米材料制備^［7］和催化化學反應^［8］等領域．

目前，大多數描寫空化泡運動的模型都是在自由場里建立的，如Rayleigh-Plesset（RP）方程^［9］和Keller-Miksis方程^［10］．考慮邊界對空化泡作用的文獻相對較少，但在實際應用中，如超聲清洗和納米材料制備總是存在不同類型的剛性壁面，因此，研究剛性壁面附近空化泡的動力學有著重要意義．鏡像法是解決點源邊值問題的有效方法之一，其基本思想是用另一個點源（像源）產生的場等效代替邊界效應．根據解存在的唯一性定理^［11］，點源和像源共同產生的場等效于需要解決的場邊值問題．本文應用鏡像法建立剛性平面附近空化氣泡的動力學方程.

1 空化泡在剛性平面附近脈動的動力學模型

1.1 鏡像原理的適用性條件

空化泡在液體聲場中脈動時，可以把空化泡看成一個聲源．根據鏡像原理，氣泡受到剛性平面的反射波可以看作虛擬鏡像氣泡的輻射波．由此，剛性界面附近空化泡的振動可以看成空化泡與其鏡像耦合成的雙氣泡振動，空化泡模型如圖1所示．

鏡像原理適用的前提是唯一性定理．為此，需要對模型作如下的假設：①空化泡在脈動過程中保持球形；②只考慮氣泡徑向振動，忽略氣泡平移運動；③考慮液體不可壓縮性；④忽略氣泡內外的質量交換以及氣泡內部的化學反應．根據假設，空化泡在剛性平面附近的脈動可以近似看成半空間內球形聲源輻射問題．輻射聲場包括兩部分^［12］：一是空化泡聲源直接到達觀測點的聲波;二是從聲源出發，經邊界面反射后再到達觀測點的聲波．由于兩個波的波陣面不同，嚴格求解這樣的聲場比較困難．目前的方法是先假設一個解，使其滿足兩個條件：一是滿足聲波動方程;二是滿足剛性表面上聲波法向速度恒為零，即紐曼邊界條件．由此，根據波動方程解的唯一性，就能確定假設解就是研究問題的唯一解．

具體過程是假設在剛性界面（分界面）的另一側，與空化泡對稱的位置上，存在振動狀況完全與原空化泡相同的另外一個聲源，用這個假定聲源輻射的聲波等效代替界面的反射聲波，把邊界問題等效為由兩個空化泡脈動源的輻射問題，其產生的合成聲壓滿足上面提到的兩個條件，即解具有惟一性．

1.2 剛性平面附近空化泡的動力學方程

用鏡像空化泡替代剛性平面后，形成的雙氣泡耦合系統的動力學可以用自由場中兩個氣泡的動力學模型來描寫^［2］，即：

在圖1中，O₁和O₂分別為空化泡與其鏡像空化泡的球心坐標，D為空化泡球心到剛性平面之間的距離．則由（1）式可以得到超聲場中剛性平面附近空化泡的動力學方程為：

在式（2）-式（5）中：p₀表示液體中的靜態壓力，一般地，p₀=p（∞）；η表示液體的黏度；σ表示空化泡的表面張力；γ表示空化泡里氣體的絕熱系數；h是表示氣泡的范德爾斯半徑，h=R₀/8.5；p_a表示驅動聲壓幅度；w=2πf，f表示聲波頻率．

2 Runge-Kutta法數值求解空化泡動力學方程

2.1 Runge-Kutta法

Runge-Kutta法是德國數學家C. Runge和M. Kutta提出的數值算法，主要用來數值求解一階常微分方程的初始值問題^［3］. 其基本思想是：從點t_i，u_i計算t_i+1，u_i+1，在計算區間t_i，t_i+1上取多個點的斜率值，然后進行加權平均作為平均斜率.考慮一個初始值問題^［4］：

其中，u₀為常數.方程（6）的泛定形式為：

在區間t，t+h內（h表示時間間隔），方程（7）可以近似表示為：

如果在區間t，t+h內取m點，即t=t₁≤t₂≤…≤t_m≤t+h.則方程（8）近似為：

其中，c_i是加權系數，滿足條件：c_i≥0且∑mi=1c_i=1.為了得到算符公式，令u（t+h）=u_n+1，u（t）=u_n.當m=1，c₁=1時，方程（9）為：

方程（10）其實就是Euler法.類似地，當m=4，c₁=16，c₂=13，c₃=13和c₄=16時，可以得到常用的四階Runge-Kutta法^［3-4］為：

2.2 空化泡動力學方程的數值解

與常規的數值求解方程步驟一樣，用Runge-Kutta數值求解空化泡的動力學方程，包括時間離散化、微分方程離散化、求解遞進方程、誤差控制和數據輸出5個過程^［3］.基于Runge-Kutta法的基本思想和算法，美國MathWorks公司在其MATrix LABorator（MATLAB）軟件包中開發了ode23和 ode45兩個常微分方程數值求解器^［13］，這里的ode是“Ordinary Different Equation”的縮寫.23和45分別表示二/三階和四/五階Runge-Kutta算法.一般地，ode45的計算精確度比ode23高，但對適度的剛性問題，ode23比ode45的計算效果好.兩者在MATLAB軟件中的調用格式均相同，即

［t，y］=ode45（‘Fun’，Tspan，y0，options）.

其中，［t，y］表示微分方程（組）的輸出數據，t和y均為列矩陣，對應每一個離散節點上的數值.Fun表示函數文件名，是輸入參數，用來定義微分方程（組）.m表示函數文件名， 引用時除了使用英文格式的單引號外，也可以在文件名前加@界定文件名.Tspan是輸入參數，表示常微分方程自變量的取值范圍，可以設置計算步長的大小，不足之處是步長設定后，在整個計算過程中是固定的，可能會造成計算時間過長.y0是輸入參數，表示初始條件向量，輸入格式為：y0=［y（t0）y′（t0）y″（t0）…］.options是輸入參數，表示選項參數，可以忽略.

為了運用Runge-Kutta法求解方程（2），需要經過以下過程：

（1）把方程（2）變成兩個一階常微分方程，即：

（2）建立以“Bubble”為文件名的子文件“Bubble.m”，子文件中包含兩個微分方程（12）和（13）；在同一個文件夾下，建立以BubDynamics為文件名的主文件“BubDynamics.m”，在主文件中調用ode45求解器數值求解子文件中的方程組：［t，y］=ode45（@Bubble，tspan，y0），其中tspan=0：step：period，period表示計算的總時間為一個聲周期，step表示計算步長，通過改變步長的大小控制計算誤差.計算中使用的物理量參數值如表1所列.

（3）使用函數save保存數據.格式為： save 'data.txt' H -ascii ，即把H矩陣中的數據保存在data.txt文本中， H=［t， y（：，1）， y（：，2）］， t表示時間列，y（：，1）表示方程（12）的數值解，即空化泡半徑隨時間的演化數據，y（：，2）表示方程（13）的數值解，即空化泡壁速度隨時間的演化數據.

（4）把data.txt 導入數據處理軟件origin8.5中進行處理，結果如圖2所示.從圖2中可以看出，當計算時間步長大于0.04 μs時，空化泡從最大值坍塌到最小值階段的相對誤差較大；當計算時間步長不大于0.04 μs時，計算結果幾乎一致，因此，本文的計算時間步長設置為0.04 μs.

3 剛性平面附近空化泡的動力學行為分析

針對文獻［2］分析了自由界面和剛性界面附近空化泡的振動，以及分析了氣泡初始半徑、氣泡到界面的距離、聲壓幅度和超聲頻率對空化泡潰滅速度的影響，但從式（2）-式（5）可以看出，除了上面4個因素影響空化泡動力學外，還有液體黏度系數和氣體絕熱系數等參數對空化泡動力學有影響．下面主要討論液體黏度系數和氣體絕熱系數分別對空化泡的振動和潰滅速度的關系．

3.1 氣體絕熱系數和液體粘度系數對空化泡徑向振動的影響

水中空氣泡、水中氬氣泡和硫酸中空氣泡在超聲作用下，其徑向半徑的時間演化如圖3所示．在剛性界面附近同一位置，在相同的超聲波驅動下，水中的空氣泡徑向振動幅度最大，氬氣泡振動幅度次之，硫酸中空氣泡振動幅度最小．這說明氣體絕熱系數越大，粘度系數越大，氣泡振動受到的抑制就越大．

征

3.2 氣體絕熱系數和液體粘度系數對空化泡潰

滅速度的影響

空化泡在材料表面塌縮時產生的微射流是超聲清洗和納米材料制備的主要機制.空化泡壁的潰滅速度與微射流密切相關．為了討論氣體絕熱系數和液體粘滯系數與空化泡潰滅速度的影響，文中分別計算了水中和硫酸中的自由空化泡，以及水-剛性平面和硫酸-剛性平面附近空化泡的潰滅速度，計算公式為^［2］：

其中：v_collapse表示空化泡壁的潰滅速度，R_max表示空化泡的最大半徑，p_v表示液體飽和蒸汽壓．

根據式（2）-式（5）和式（14），計算得到水中和硫酸中的自由空化泡、水-剛性平面和硫酸-剛性平面附近空化泡的潰滅速度，其關系如圖4所示.

從圖4可以看出不論是自由場中的空化泡，還是剛性平面附近的空化泡，不論在水中還是硫酸中，空化泡潰滅速度均隨著氣體絕熱系數的增大而減小；硫酸中空化泡的潰滅速度遠小于水中的空化泡的潰滅速度．這說明液體粘度系數越大（硫酸的粘度系數比水的大），空化泡的潰滅速度越小．

4 結論

本文首先用鏡像法得到了剛性平面附近空化氣泡的動力學方程，然后基于Runge-Kutta法，用MATLAB軟件包編程計算了該方程，分析了氣體絕熱系數和液體粘度系數對空化泡徑向振動及其潰滅速度的影響．結果表明氣體絕熱系數和粘度系數越大，空化泡振動受到的抑制就越大，空化泡的潰滅速度越小．相同的液體環境和超聲驅動下，自由場中空化泡的潰滅速度均比邊界場中空化泡的潰滅速度大．這些計算結果不僅表明了四階Runge-Kutta法在求解非線性常微分方程的可靠性，而且有利于理解空化氣泡與界面的相互作用機制，對超聲清洗、超聲制備材料等方面的應用具有指導意義.

參考文獻：

［1］ LIANG J F，CHEN W Z，SHAO W H，et al.Aspherical oscillation of two interacting bubbles in an ultrasound field ［J］.Chinese Physics Letters，2012，29（7）：074701.

［2］ 郭策，祝錫晶，王建青，等.超聲場下剛性界面附近潰滅空化氣泡的速度分析［J］.物理學報，2016，65（4）：184-191.

［3］ 陳偉中.聲空化物理［M］.北京：高等教育出版社，2014.

［4］ 李榮華，劉播.微分方程數值解法 ［M］.第4版.北京：高等教育出版社，2019.

［5］ 安宇.空化氣泡發光機理［J］.中國科學：物理學力學，2011，41（4）：343-349．

［6］ NIEMCZEWSKI B.Cavitation intensity of water under practical ultrasonic cleaning conditions［J］.Ultrasonics Sonochemistry，2014，21：354-359.

［7］ OSHIMA R，YAMAMOTO T A，MIZUKOSHI Y，et al.Electronmicroscopy of noble metal alloy nanoparticles prepared by sonochemical methods［J］.Nanostructured Materials，1999，12：111-114.

［8］ HUSSEIN E M，KHAIROU K S.Sonochemistry：Synthesis of bioactive heterocycles［J］.Review Journal of Chemistry，2014，4：221-251.

［9］ PLESSET M S.On the stability of fluid flows with spherical symmetry［J］.Journal of Applied Physics，1954，25（1）：96-98.

［10］ KELLER J B，MIKSIS M.Bubble oscillations of large amplitude［J］.The Journal of the Acoustical Society of America，1980，68（2）：628-633.

［11］ 郭碩鴻.電動力學 ［M］.第3版.北京：高等教育出版社，2008.

［12］ 杜功煥，朱哲民，龔秀芬.聲學基礎 ［M］.第2版.南京：南京大學出版社，2001.

［13］ 陸君安，尚濤.偏微分方程的MATLAB解法［M］.武漢：武漢大學出版社，2001.

［責任編輯：李 嵐］