基于Matlab GUI的气泡动力学仿真系统设计

王 寻， 王宏哲， 张泽坤， 周 敏， 梁金福

（1.上海电机学院a.凯撒斯劳滕智能制造学院；b.文理学院，上海 201306；2.西安工程大学理学院，西安 710048；3.贵州师范大学物理与电子科学学院，贵阳 550025）

0 引 言

当大功率超声在液体中传播时，会使得液体中原先存在的气核逐渐增大为肉眼可见的气泡。这种现象称为声空化［1］。声空化在催化［2］、清洗［3］和乳化［4］等领域都有广泛的应用，这些应用主要得益于空化泡效应。空化泡在超声作用下的运动过程较为复杂，在平移的同时也会发生高频振动［5］。若能将气泡运动以三维动态图像的形式展现出来，则有利于学生直观了解气泡在超声驱动下的动力学行为，更好地掌握超声空化相关知识。

Matlab是目前数值计算和数据分析领域的主流软件之一，具有丰富的矩阵运算和绘图功能。Matlab GUI具有功能强大、简单易学等优点，能便捷地实现用户交互，将计算结果以图形化的形式展示给用户，在科学研究和教学中得到了广泛的应用。已有很多学者基于Matlab GUI设计了仿真软件，例如张义灵等［6］使用Matlab软件对水平面上受不稳定约束的弹簧振子动力学方程进行数值求解，并通过Matlab GUI设计了图形界面。用户可在图形界面中输入弹簧劲度系数、弹簧原始长度和振子质量等参数，观察弹簧振子随时间变化的曲线以及弹簧振子的运动轨迹。陈梅等［7］基于Matlab GUI 实现了Ziegler-Nichols PID 参数整定的仿真系统。该系统能够显示单位阶跃响应曲线，便于用户快速直观地获取PID参数。孙少华等［8］基于Matlab GUI设计了交互式混合动力汽车教学仿真平台。利用该平台，用户可以实现混合动力驱动系统方案选择、动力传动系统匹配和整车及核心部件参数设置等功能。该平台界面友好，操作便捷，可有效促进学生对理论知识的掌握。本研究基于Matlab GUI 设计了图形界面软件，软件将根据用户输入的驱动声压幅值、驱动声波频率、液体黏滞系数和气泡平衡半径等参数，使用Matlab对气泡动力学方程进行数值求解，并通过三维图像直观展示气泡的脉动过程。该软件可通过交互过程激发学生的学习兴趣，改善教学效果，增强学生对相关内容的理解。

1 气泡动力学的理论基础

常用于描述超声作用下气泡动力学的模型包括Rayleigh-Plesset（R-P）模型［9-10］、Keller-Miksis（K-M）模型［11-13］和Gilmore 模型［14］等。K-M 模型考虑了液体的可压缩性，能够更加真实地反映气泡动力学规律，本文所述的软件开发中采用K-M模型。

从气泡数量上看，气泡动力学问题主要可被归类为单气泡动力学、双气泡动力学和气泡群动力学。本文将实现单气泡与双气泡的数值计算及气泡二维和三维动态图像演示。

1.1 单气泡动力学

在考虑平移的情况下，单个气泡在声波作用下的动力学方程组［15］：

式中：R为气泡半径；x为气泡位置；c为液体中声速；ρ为液体密度；Psc为气泡表面液体的压强；Fex为作用在气泡上的外力。

需要指出的是，本文仅考虑气泡沿x轴运动的情况。Fex可以表示为Bjerknes 力（Fpr）和黏滞力（Fvis）之和，即

Fpr和Fvis可分别表示为：

式中：Pa、ω 分别为驱动声波的声压幅度和角频率；k为波数；η为液体黏滞系数；d为气泡中心的x坐标与距离最近的声压波腹的x坐标之差，设定与气泡中心距离最近的声压波腹x=0。

式（1）中Psc可以表示为

式中：R0为气泡的平衡半径；P0为液体的静态压强；σ为表面张力系数；γ 为气体绝热系数；Pex为驱动声压［11］，

联立式（1）～（6），可通过数值计算得出单气泡半径和位置随时间的变化。

1.2 双气泡动力学

当声波驱动2 个气泡组成的系统时，气泡间存在次Bjerknes力的相互作用，任意一个气泡的平移和脉动都会受到另一个气泡的影响。假定2 个气泡都在声压波腹点附近，气泡动力学方程组为［16］：

式中：Ri为第i个气泡的半径；D为两气泡间距。两个气泡表面的压强可以表示为

式中：Ri0为第i个气泡的平衡半径；Pexi为作用在第i个气泡上的驱动声压。

气泡受到的外力

式中，υ3-i为第3 -i个气泡在第i个气泡中心位置产生的液体流速，可表示为

在本研究中，设置驱动声压为正弦波

联立式（7）～（11），可通过数值计算得出两个气泡半径和位置随时间的变化。

2 气泡动力学方程的数值求解

为便于数值计算，需对气泡动力学方程组进行降阶处理［6］。当求解单气泡动力学问题时，R和x为待求解的变量。结合式（1）～（6）可以表示出d2R/dt2和d2x/dt2，构建向量y =［dR/dt，d2R/dt2，dx/dt，d2x/dt2］。将y代入Matlab中ode45 函数，基于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。得到的结果如图1 所示。由图1（a）可见，气泡半径呈周期性变化，在每一个脉动周期中，气泡都会经历膨胀和收缩的过程；图1（b）表明，气泡将在1/4 波长位置（x=0.25λ）附近抖动；结合图1（a）、（b）可见，气泡脉动时半径变化的幅值与气泡位置有关，随着气泡位置的改变，气泡半径变化幅值也呈周期性变化。

图1 单气泡的半径和位置随时间的变化

求解双气泡问题时，待求解变量为R1、R2、x1、x2。结合式（7）～（11），可构建向量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 函数进行数值计算，得到2 个气泡半径和位置变化曲线如图2（a）、（b）所示。由图2（b）可见，两气泡先相互靠近，然后保持固定的距离移动；由图2（a）可见，在整个仿真过程中，气泡1脉动时半径变化幅值几乎保持不变。在相互靠近阶段，气泡2 半径幅值逐渐减小，在保持固定距离移动阶段，气泡2 半径幅值保持不变。这与文献［16］中的研究是一致的。

图2 双气泡的半径和位置随时间的变化

3 基于Matlab GUI的仿真界面设计

3.1 总体设计

为实现用户与计算机的便捷交互，在实现气泡动力学数值计算程序后，需进行图形界面设计。本研究设计的软件主界面如图3 所示。通过选择“单气泡演示”或“双气泡演示”按钮，可跳转至相应界面，以便于用户进行仿真。

3.2 单气泡演示图形界面

如图4 所示，单气泡界面可分为参数输入区域和图形显示区域。参数输入区域中的参数包括声速、液体密度、驱动声压幅度、仿真时间等。图形显示区域分别展示气泡半径和位置变化曲线以及气泡运动的动态过程。图形可放大或缩小。

图3 气泡动力学仿真软件主界面

图4 单气泡运动演示界面

用户在文本框中输入的数据，需要通过程序读取并用于数值计算。例如驱动频率，由对应的edit 控件名称为edit_fre，软件中需使用

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

语句将用户在控件中输入的值读入‘fre’变量中。其他变量读取与此类似。在读取完成后，调用ode45 函数进行数值求解。求解得到的半径和位置变量存放在‘r1’和‘p1’变量中。通过plot 函数画出‘r1’和‘p1’随时间的变化（图4 中右上部分）。

值得重点介绍的是三维动态图像显示的实现。由于软件中ode45 求解结果中的‘t’并非均匀分布，为了画图时能够匀速显示气泡的平移和脉动过程，需要对‘t’和其对应‘r1’和‘p1’的数据进行等间隔重采样。方法如下：

其中axis指令用于指定三维坐标轴的范围，需要使得x方向坐标轴能够包括气泡的平移范围，且必须确保气泡膨胀到最大时仍能被完整显示。‘rsout. Data’和‘psout.Data’为气泡半径和位置数据。sphere（80）表示产生一个由80 ×80 个面组成的球体，surf为球体绘制函数。colormap指令控制球体颜色映射，camlight控制图像光源，draw用于刷新图像，cla用于清除图像。

3.3 双气泡演示图形界面

在图3 界面中点击“双气泡演示”时，跳转到双气泡动力学演示界面（见图5）。界面与图4 所示的单气泡动力学演示界面相似，但参数有所增加。其中2 个气泡的平衡半径和初始速度都有待于用户输入。另外还多了一个按钮组，让用户选择显示二维图像还是三维图像。此界面中，所有展示的图像都包含了2 个气泡的信息。

图5 双气泡运动展示界面

由于按钮组中的单选按钮同时只能有一个处于选中状态，软件中通过

指令读取“二维演示”这个单选按钮当前是否选中，并通过

控制显示的动态图为二维还是三维。若“二维演示”未被选中，则调用view（3）函数。若被选中则调用view（2）函数。图5 中右下部分为气泡脉动三维演示。

当实现双气泡动力学动态演示时，由于要同时显示2 个气泡的脉动和平移，须对单气泡演示的程序进行改写。要对计算得到的2 个气泡半径和位置数据进行重采样：

其中hold on语句可以让第1 个球体画完后保持显示，再画第2 个球体。其余函数含义与3.2 节所述单气泡图像绘制相同。

若使用二维演示，则显示界面如图6 所示。此界面下方显示的为气泡的二维动态图像。二维图像的优点在用户可以方便地读取气泡位置坐标和半径，但其直观性可能不如三维图像。

图6 双气泡动力学的二维动态图演示

4 结 语

本文基于Matlab GUI 设计了超声驱动下气泡动力学仿真系统。软件允许用户自行输入驱动声压幅值、驱动声波频率、液体黏滞系数和气泡平衡半径等参数，进行单气泡和多气泡动力学的数值计算，并将计算结果以曲线形式显示在界面中。为增强结果显示的直观效果，设计了动态图显示功能，将气泡的脉动和平移过程动态显示在界面上。该软件可通过与学生的交互，提高学生的学习兴趣，增强学生对声空化相关知识的理解。未来将组织物理专业学生进行实验，定量研究该软件在声空化教学中的有效性，并根据用户反馈对软件进行改进，不断完善软件功能。