王 尋， 王宏哲， 張澤坤， 周 敏， 梁金福

（1.上海電機學院a.凱撒斯勞滕智能制造學院；b.文理學院，上海 201306；2.西安工程大學理學院，西安 710048；3.貴州師范大學物理與電子科學學院，貴陽 550025）

0 引 言

當大功率超聲在液體中傳播時，會使得液體中原先存在的氣核逐漸增大為肉眼可見的氣泡。這種現(xiàn)象稱為聲空化［1］。聲空化在催化［2］、清洗［3］和乳化［4］等領域都有廣泛的應用，這些應用主要得益于空化泡效應。空化泡在超聲作用下的運動過程較為復雜，在平移的同時也會發(fā)生高頻振動［5］。若能將氣泡運動以三維動態(tài)圖像的形式展現(xiàn)出來，則有利于學生直觀了解氣泡在超聲驅(qū)動下的動力學行為，更好地掌握超聲空化相關知識。

Matlab是目前數(shù)值計算和數(shù)據(jù)分析領域的主流軟件之一，具有豐富的矩陣運算和繪圖功能。Matlab GUI具有功能強大、簡單易學等優(yōu)點，能便捷地實現(xiàn)用戶交互，將計算結(jié)果以圖形化的形式展示給用戶，在科學研究和教學中得到了廣泛的應用。已有很多學者基于Matlab GUI設計了仿真軟件，例如張義靈等［6］使用Matlab軟件對水平面上受不穩(wěn)定約束的彈簧振子動力學方程進行數(shù)值求解，并通過Matlab GUI設計了圖形界面。用戶可在圖形界面中輸入彈簧勁度系數(shù)、彈簧原始長度和振子質(zhì)量等參數(shù)，觀察彈簧振子隨時間變化的曲線以及彈簧振子的運動軌跡。陳梅等［7］基于Matlab GUI 實現(xiàn)了Ziegler-Nichols PID 參數(shù)整定的仿真系統(tǒng)。該系統(tǒng)能夠顯示單位階躍響應曲線，便于用戶快速直觀地獲取PID參數(shù)。孫少華等［8］基于Matlab GUI設計了交互式混合動力汽車教學仿真平臺。利用該平臺，用戶可以實現(xiàn)混合動力驅(qū)動系統(tǒng)方案選擇、動力傳動系統(tǒng)匹配和整車及核心部件參數(shù)設置等功能。該平臺界面友好，操作便捷，可有效促進學生對理論知識的掌握。本研究基于Matlab GUI 設計了圖形界面軟件，軟件將根據(jù)用戶輸入的驅(qū)動聲壓幅值、驅(qū)動聲波頻率、液體黏滯系數(shù)和氣泡平衡半徑等參數(shù)，使用Matlab對氣泡動力學方程進行數(shù)值求解，并通過三維圖像直觀展示氣泡的脈動過程。該軟件可通過交互過程激發(fā)學生的學習興趣，改善教學效果，增強學生對相關內(nèi)容的理解。

1 氣泡動力學的理論基礎

常用于描述超聲作用下氣泡動力學的模型包括Rayleigh-Plesset（R-P）模型［9-10］、Keller-Miksis（K-M）模型［11-13］和Gilmore 模型［14］等。K-M 模型考慮了液體的可壓縮性，能夠更加真實地反映氣泡動力學規(guī)律，本文所述的軟件開發(fā)中采用K-M模型。

從氣泡數(shù)量上看，氣泡動力學問題主要可被歸類為單氣泡動力學、雙氣泡動力學和氣泡群動力學。本文將實現(xiàn)單氣泡與雙氣泡的數(shù)值計算及氣泡二維和三維動態(tài)圖像演示。

1.1 單氣泡動力學

在考慮平移的情況下，單個氣泡在聲波作用下的動力學方程組［15］：

式中：R為氣泡半徑；x為氣泡位置；c為液體中聲速；ρ為液體密度；Psc為氣泡表面液體的壓強；Fex為作用在氣泡上的外力。

需要指出的是，本文僅考慮氣泡沿x軸運動的情況。Fex可以表示為Bjerknes 力（Fpr）和黏滯力（Fvis）之和，即

Fpr和Fvis可分別表示為：

式中：Pa、ω 分別為驅(qū)動聲波的聲壓幅度和角頻率；k為波數(shù)；η為液體黏滯系數(shù)；d為氣泡中心的x坐標與距離最近的聲壓波腹的x坐標之差，設定與氣泡中心距離最近的聲壓波腹x=0。

式（1）中Psc可以表示為

式中：R0為氣泡的平衡半徑；P0為液體的靜態(tài)壓強；σ為表面張力系數(shù)；γ 為氣體絕熱系數(shù)；Pex為驅(qū)動聲壓［11］，

聯(lián)立式（1）～（6），可通過數(shù)值計算得出單氣泡半徑和位置隨時間的變化。

1.2 雙氣泡動力學

當聲波驅(qū)動2 個氣泡組成的系統(tǒng)時，氣泡間存在次Bjerknes力的相互作用，任意一個氣泡的平移和脈動都會受到另一個氣泡的影響。假定2 個氣泡都在聲壓波腹點附近，氣泡動力學方程組為［16］：

式中：Ri為第i個氣泡的半徑；D為兩氣泡間距。兩個氣泡表面的壓強可以表示為

式中：Ri0為第i個氣泡的平衡半徑；Pexi為作用在第i個氣泡上的驅(qū)動聲壓。

氣泡受到的外力

式中，υ3-i為第3 -i個氣泡在第i個氣泡中心位置產(chǎn)生的液體流速，可表示為

在本研究中，設置驅(qū)動聲壓為正弦波

聯(lián)立式（7）～（11），可通過數(shù)值計算得出兩個氣泡半徑和位置隨時間的變化。

2 氣泡動力學方程的數(shù)值求解

為便于數(shù)值計算，需對氣泡動力學方程組進行降階處理［6］。當求解單氣泡動力學問題時，R和x為待求解的變量。結(jié)合式（1）～（6）可以表示出d2R/dt2和d2x/dt2，構(gòu)建向量y =［dR/dt，d2R/dt2，dx/dt，d2x/dt2］。將y代入Matlab中ode45 函數(shù)，基于4 階5級Runge-Kutta變步長算法求解氣泡動力學方程組，即可得到氣泡脈動過程中半徑和位置隨時間的變化。求解時，設置：R0=93 μm，f=25 kHz，ρ =998 kg/m3，c=1 500 m/s，P0=0.101 3 MPa，Pa=0.13 MPa，γ=1.4，σ=72.5 mN/m，μ=1 mPa·s。得到的結(jié)果如圖1 所示。由圖1（a）可見，氣泡半徑呈周期性變化，在每一個脈動周期中，氣泡都會經(jīng)歷膨脹和收縮的過程；圖1（b）表明，氣泡將在1/4 波長位置（x=0.25λ）附近抖動；結(jié)合圖1（a）、（b）可見，氣泡脈動時半徑變化的幅值與氣泡位置有關，隨著氣泡位置的改變，氣泡半徑變化幅值也呈周期性變化。

圖1 單氣泡的半徑和位置隨時間的變化

求解雙氣泡問題時，待求解變量為R1、R2、x1、x2。結(jié)合式（7）～（11），可構(gòu)建向量y =［dR1/dt，d2R1/dt2，dR2/dt，d2R2/dt2，dx1/dt，d2x1/dt2，dx2/dt，d2x2/dt2］。設置：R10=4 μm，R20=3 μm，μ =1.0 mPa·s，f=20 kHz，L=0.2 mm，σ=72.5 mN/m，P0=0.101 3 MPa，Pa=1.21 P。將y 代入ode45 函數(shù)進行數(shù)值計算，得到2 個氣泡半徑和位置變化曲線如圖2（a）、（b）所示。由圖2（b）可見，兩氣泡先相互靠近，然后保持固定的距離移動；由圖2（a）可見，在整個仿真過程中，氣泡1脈動時半徑變化幅值幾乎保持不變。在相互靠近階段，氣泡2 半徑幅值逐漸減小，在保持固定距離移動階段，氣泡2 半徑幅值保持不變。這與文獻［16］中的研究是一致的。

圖2 雙氣泡的半徑和位置隨時間的變化

3 基于Matlab GUI的仿真界面設計

3.1 總體設計

為實現(xiàn)用戶與計算機的便捷交互，在實現(xiàn)氣泡動力學數(shù)值計算程序后，需進行圖形界面設計。本研究設計的軟件主界面如圖3 所示。通過選擇“單氣泡演示”或“雙氣泡演示”按鈕，可跳轉(zhuǎn)至相應界面，以便于用戶進行仿真。

3.2 單氣泡演示圖形界面

如圖4 所示，單氣泡界面可分為參數(shù)輸入?yún)^(qū)域和圖形顯示區(qū)域。參數(shù)輸入?yún)^(qū)域中的參數(shù)包括聲速、液體密度、驅(qū)動聲壓幅度、仿真時間等。圖形顯示區(qū)域分別展示氣泡半徑和位置變化曲線以及氣泡運動的動態(tài)過程。圖形可放大或縮小。

圖3 氣泡動力學仿真軟件主界面

圖4 單氣泡運動演示界面

用戶在文本框中輸入的數(shù)據(jù)，需要通過程序讀取并用于數(shù)值計算。例如驅(qū)動頻率，由對應的edit 控件名稱為edit_fre，軟件中需使用

fre = str2num（get（handles.edit_fre，'String'））；

語句將用戶在控件中輸入的值讀入‘fre’變量中。其他變量讀取與此類似。在讀取完成后，調(diào)用ode45 函數(shù)進行數(shù)值求解。求解得到的半徑和位置變量存放在‘r1’和‘p1’變量中。通過plot 函數(shù)畫出‘r1’和‘p1’隨時間的變化（圖4 中右上部分）。

值得重點介紹的是三維動態(tài)圖像顯示的實現(xiàn)。由于軟件中ode45 求解結(jié)果中的‘t’并非均勻分布，為了畫圖時能夠勻速顯示氣泡的平移和脈動過程，需要對‘t’和其對應‘r1’和‘p1’的數(shù)據(jù)進行等間隔重采樣。方法如下：

其中axis指令用于指定三維坐標軸的范圍，需要使得x方向坐標軸能夠包括氣泡的平移范圍，且必須確保氣泡膨脹到最大時仍能被完整顯示。‘rsout. Data’和‘psout.Data’為氣泡半徑和位置數(shù)據(jù)。sphere（80）表示產(chǎn)生一個由80 ×80 個面組成的球體，surf為球體繪制函數(shù)。colormap指令控制球體顏色映射，camlight控制圖像光源，draw用于刷新圖像，cla用于清除圖像。

3.3 雙氣泡演示圖形界面

在圖3 界面中點擊“雙氣泡演示”時，跳轉(zhuǎn)到雙氣泡動力學演示界面（見圖5）。界面與圖4 所示的單氣泡動力學演示界面相似，但參數(shù)有所增加。其中2 個氣泡的平衡半徑和初始速度都有待于用戶輸入。另外還多了一個按鈕組，讓用戶選擇顯示二維圖像還是三維圖像。此界面中，所有展示的圖像都包含了2 個氣泡的信息。

圖5 雙氣泡運動展示界面

由于按鈕組中的單選按鈕同時只能有一個處于選中狀態(tài)，軟件中通過

指令讀取“二維演示”這個單選按鈕當前是否選中，并通過

控制顯示的動態(tài)圖為二維還是三維。若“二維演示”未被選中，則調(diào)用view（3）函數(shù)。若被選中則調(diào)用view（2）函數(shù)。圖5 中右下部分為氣泡脈動三維演示。

當實現(xiàn)雙氣泡動力學動態(tài)演示時，由于要同時顯示2 個氣泡的脈動和平移，須對單氣泡演示的程序進行改寫。要對計算得到的2 個氣泡半徑和位置數(shù)據(jù)進行重采樣：

其中hold on語句可以讓第1 個球體畫完后保持顯示，再畫第2 個球體。其余函數(shù)含義與3.2 節(jié)所述單氣泡圖像繪制相同。

若使用二維演示，則顯示界面如圖6 所示。此界面下方顯示的為氣泡的二維動態(tài)圖像。二維圖像的優(yōu)點在用戶可以方便地讀取氣泡位置坐標和半徑，但其直觀性可能不如三維圖像。

圖6 雙氣泡動力學的二維動態(tài)圖演示

4 結(jié) 語

本文基于Matlab GUI 設計了超聲驅(qū)動下氣泡動力學仿真系統(tǒng)。軟件允許用戶自行輸入驅(qū)動聲壓幅值、驅(qū)動聲波頻率、液體黏滯系數(shù)和氣泡平衡半徑等參數(shù)，進行單氣泡和多氣泡動力學的數(shù)值計算，并將計算結(jié)果以曲線形式顯示在界面中。為增強結(jié)果顯示的直觀效果，設計了動態(tài)圖顯示功能，將氣泡的脈動和平移過程動態(tài)顯示在界面上。該軟件可通過與學生的交互，提高學生的學習興趣，增強學生對聲空化相關知識的理解。未來將組織物理專業(yè)學生進行實驗，定量研究該軟件在聲空化教學中的有效性，并根據(jù)用戶反饋對軟件進行改進，不斷完善軟件功能。