基于Matlab GUI的氣泡動(dòng)力學(xué)仿真系統(tǒng)設(shè)計(jì)

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

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

0 引 言

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

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

1 氣泡動(dòng)力學(xué)的理論基礎(chǔ)

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

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

1.1 單氣泡動(dòng)力學(xué)

在考慮平移的情況下，單個(gè)氣泡在聲波作用下的動(dòng)力學(xué)方程組［15］：

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

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

Fpr和Fvis可分別表示為：

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

式（1）中Psc可以表示為

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

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

1.2 雙氣泡動(dòng)力學(xué)

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

式中：Ri為第i個(gè)氣泡的半徑；D為兩氣泡間距。兩個(gè)氣泡表面的壓強(qiáng)可以表示為

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

氣泡受到的外力

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

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

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

2 氣泡動(dòng)力學(xué)方程的數(shù)值求解

為便于數(shù)值計(jì)算，需對(duì)氣泡動(dòng)力學(xué)方程組進(jìn)行降階處理［6］。當(dāng)求解單氣泡動(dòng)力學(xué)問題時(shí)，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級(jí)Runge-Kutta變步長算法求解氣泡動(dòng)力學(xué)方程組，即可得到氣泡脈動(dòng)過程中半徑和位置隨時(shí)間的變化。求解時(shí)，設(shè)置：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）可見，氣泡半徑呈周期性變化，在每一個(gè)脈動(dòng)周期中，氣泡都會(huì)經(jīng)歷膨脹和收縮的過程；圖1（b）表明，氣泡將在1/4 波長位置（x=0.25λ）附近抖動(dòng)；結(jié)合圖1（a）、（b）可見，氣泡脈動(dòng)時(shí)半徑變化的幅值與氣泡位置有關(guān)，隨著氣泡位置的改變，氣泡半徑變化幅值也呈周期性變化。

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

求解雙氣泡問題時(shí)，待求解變量為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］。設(shè)置：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ù)進(jìn)行數(shù)值計(jì)算，得到2 個(gè)氣泡半徑和位置變化曲線如圖2（a）、（b）所示。由圖2（b）可見，兩氣泡先相互靠近，然后保持固定的距離移動(dòng)；由圖2（a）可見，在整個(gè)仿真過程中，氣泡1脈動(dòng)時(shí)半徑變化幅值幾乎保持不變。在相互靠近階段，氣泡2 半徑幅值逐漸減小，在保持固定距離移動(dòng)階段，氣泡2 半徑幅值保持不變。這與文獻(xiàn)［16］中的研究是一致的。

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

3 基于Matlab GUI的仿真界面設(shè)計(jì)

3.1 總體設(shè)計(jì)

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

3.2 單氣泡演示圖形界面

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

圖3 氣泡動(dòng)力學(xué)仿真軟件主界面

圖4 單氣泡運(yùn)動(dòng)演示界面

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

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

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

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

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

3.3 雙氣泡演示圖形界面

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

圖5 雙氣泡運(yùn)動(dòng)展示界面

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

指令讀取“二維演示”這個(gè)單選按鈕當(dāng)前是否選中，并通過

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

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

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

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

圖6 雙氣泡動(dòng)力學(xué)的二維動(dòng)態(tài)圖演示

4 結(jié) 語

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