空化水射流空泡溃灭过程的数值分析

杨博凯，卢义玉，杨晓峰，冯明涛，张 赛

(重庆大学资源及环境科学学院，重庆400030)

0 引言

目前，高压空化水射流已经应用于很多领域，如岩石破碎、海上石油钻井、船舶清洗、空化水射流处理有机物废水以及医疗技术等，但目前对于空泡生成和溃灭的机理还不清楚.

对于空化水射流而言，其空泡半径与压力、时间存在确定的函数关系，而三者的关系直接影响空化射流的效率以及能量释放程度.因此，旨在研究它们之间关系的空泡动力学方程在空化理论中占有重要的地位.

对于一个包含空泡的液体，如果已知射流中单个空泡的内外压力，就可以计算它的半径R［1］，同时能够建立一个Rayleigh-Plesset方程来近似求得气泡的运动状态［2］.理论上这个方程可以使用常见的任意一种有限元差分法解决，但由于空泡半径随时间变化的过程是非线性的，常用的迭代方法在求解Rayleigh-Plesset方程时存在困难，计算量大、精确度低，特别是在空泡溃灭阶段会产生振荡.笔者根据变步长法的基本思想，结合自主研发的空化水射流对各种迭代法的计算效率进行分析，提出了变步长的Rayleigh-Plesset方程计算方法.变步长法能够更有效地模拟空泡最终溃灭阶段的能量释放过程，而这个能量释放过程正是研究其对附近边界产生破坏作用机理的关键［3－4］.

1 Rayleigh-Plesset空泡运动方程

Rayleigh-Plesset方程是解决空泡运动方程的经典方程，在空化水射流的研究中，经常需要确定空泡随时间的变化以及最终溃灭放能情况，而通过Rayleigh-Plesset方程可以较为准确地得到空泡半径与时间的关系.它是一个二阶微分方程［5］.

式中：ρ为液体密度;μ为液体运动黏度;t为时间;p∞为液体静压;pB为空泡内的压力.

式(1)可以转换为以下的一阶方程

当x为因变量时，令R=x，可同时求解式(2)和(3)以得到方程的解.令：

采用 Euler法、Central法、改进 Euler法和Runge－Kutta－Fehlberg法4种有限元拆分法求解［6－7］，则需将方程(1)改写成

先将式(7)转变成形如式(2)和式(3)的两个一阶等式，以便应用Runge－Kutta－Fehlberg法，令 x=R，则

假定 pB和 p∞以及初始条件已知，便可以通过式(8)和(9)求解空泡半径R，即空泡从生长到溃灭的过程中的大小可以以时间的函数来表示.

在式(9)中，变量是pB和p∞.推测空泡中含有杂质气体，其压力为pgo大小为Ro，水蒸气压力为pv，在气体不会凝结的前提下，它是一个与温度相关的定值［9－10］，因此，笔者将其简化为

则空泡压力可以表示为

液体压力p∞(t)是一个重要的参数，它关系到空泡成长和溃灭的关键参数p∞.没有可以直接求得p∞的准确表达式，一般由经验公式来获得，p∞由CFD软件模拟的结果获得.

为解Rayleigh-Plesset方程，需要先确定某些参数.笔者的测试条件为产生空泡的介质为水，且温度保持在300 K不变，水的性能见表1.

表1 300 K时水的性能Tab.1 Water properties at 300 K

2 固定步长法解Rayleigh-Plesset方程

2.1 固定步长法的基本算法

在满足前述解Rayleigh-Plesset方程各项假设的前提下，对于固定步长法的计算方法可表示如下：(1)令 h0=Δt;(2)取 ti=ti－1+h0i;(3)将结果带入式(8)和(9)中;(4)由Rayleigh-Plesset方程求得时刻的空泡半径.

2.2 几种数值解法的稳定性分析

Rayleigh-Plesset方程是一个非线性微分方程，且没有解析解，不易比较判断各种数值方法的可靠性，因此采用下式来进行测试.

式(12)在x=90°时具有奇异性，与原式性质相同，因此符合要求［8］.应用式(12)并利用 EXCEL工具，得到的4种数值方法从Δx=1°开始的步长算法分析如图1.

图1 对比分析几种方程解的情况Fig.1 Comparison of solutions from different schemes with analytical solution

如图1所示，在x不接近90°时，各种方法都能得到较为准确的结果，但在整个0°≤x≤90°的范围中考虑的话，Runge-Kutta-Fehlberg法效果最好.

4种方法中只有Runge-Kutta-Fehlberg法可以估算误差.对比解析解可以看出Runge-Kutta-Fehlberg法的估算误差较为保守一些.

2.3 较小压力下各种解法解Rayleigh-Plesset方程

如表2所示，压力从12 kPa线性下降到最低值－1 kPa，之后又恢复到12 kPa.时间步长h=Δt=2 ×10－9s，18 500 步后，4 种方法的解如图2.

表2 较小压力变化下的测试情况Tab.2 Assumed pressure for testing numerical methods

图2 在固定时间步长为2×10－9s时不同方法求解Rayleigh-Plesset方程的情况Fig 2 Solutions for Rayleigh-Plesset equation by different methods with constant time step of second

从图2中可以看出，除Central法在空泡溃灭再循环的时候产生严重的数值振荡之外，其他三种方法得到了稳定且几乎是相同的解.

2.4 压力变化较大时解Rayleigh-Plesset方程

根据K.BREMHORST等人的研究可知，在较大压力变化的情况下，所需的时间步长非常小.因此，扩大压力的变化范围进行实验，使最高和最低压力分别为120 kPa和－10 kPa，如表3所示.当时间步长为2×10－9s时，即使在一个生长溃灭循环之内，4种方法也都不能完成.

表3 扩大压力变化时的测试情况Tab.3 Enlarged pressure variation for testing numerical methods

为尽量完成至少一个循环，缩减时间步长至1×20－9s再次计算，但从结果上看不出有明显的变化.这是因为当空泡即将溃灭时，时间步长与空泡半径达到最小值的变化率相比过大，所以除非大幅缩减时间步长，否则不会有明显的改观.

4种方法在溃灭结束、开始第二次循环的时候都出现了振荡现象，没有给出合适的解.这表明，由于溃灭的最后阶段压力变化的速度和空泡半径的变化速度都很快，固定时间步长的方法解Rayleigh-Plesset方程需要一个极小的时间步长，而过小的步长将导致计算步数与计算时间成几何倍数增长，这种增长往往超过了现有计算机所能承受的极限.因此，必须找到一种更为合理的算法，才能在压力变化比较大的情况下较好的解决Rayleigh-Plesset方程.

3 变步长法解Rayleigh-Plesset方程

3.1 变步长法的基本算法

在满足前述解Rayleigh-Plesset方程各项假设的前提下，对于变步长法的计算方法可表示如下：(1)令 h0= Δt，hi=hi－1Ri/Ri－1;(2)取 ti=ti－1+hi;(3)将结果带入式(8)和式(9)中;(4)由Rayleigh-Plesset方程求得ti时刻的空泡半径Ri.

3.2 可变步长法及其在压力变化较大时的测试

Rayleigh-Plesset方程的特点是可以描述空泡生长至溃灭的全过程，当空泡半径很小时，它的导数变得很大，因此需要非常小的时间步长.从上述计算情况来看，对于二阶Rayleigh-Plesset微分方程，选取的时间步长并不合适.

基于上述考虑，提出了一种时间步长随空泡半径变化速率Ri/Ri－1而自动调整的算法，配合一些简单的算法，如Euler法，来达到变步长解Rayleigh-Plesset方程的目的.图5显示了使用这种变步长法在压力变化较大时测得的结果，与从图2到图4的结果相比，这无疑是一个更为准确的结果.更重要的是，它在的时候并不会产生奇异点，因此可以连续地表示一个循环结束到下一循环开始的情况.使用这种变步长结合Euler法，计算空泡4次生长溃灭循环仅需不到1 500步.图6是空泡半径R随时间变化的情况.图中显示，空泡结束第一次循环，完成溃灭并开始反弹的时间是t=0.000 015 s，这与固定步长法测得的结果相同，但固定步长法无法模拟溃灭结束并开始反弹的过程.

3.3 变步长法与固定步长法的比较

为对比变步长法与固定步长法的计算效果，使用固定步长法中较精确的Runge-Kutta-Fehlberg法与变步长法进行对比.

在较小压力变化的情况下，分别采用变步长配合Euler法与Runge-Kutta-Fehlberg法对一个生长溃灭循环计算，其中 pmax=12 kPa，pmin=－1 kPa.结果见图 7.

如图所示，两种方法的结果没有显著的差异，均可模拟在压力变化不大情况下溃灭和反弹的过程.然而在迭代步数上却有很大差异.Runge-Kutta-Fehlberg法共使用16 000步，而变步长配合Euler法仅使用800步就得到了正确的结果.

图7 固定Runge-Kutta法与变步长配合Euler法在较低压力变化下的比较Fig.7 Comparison of solutions from the constant time step Runge-Kutta method and the variable time step on Euler method under the smaller pressure variation

在变步长法中，最大和最小的时间增量分别为10－7s和10－17s.若使用固定步长法，则必须采用最小的步长以达到相同的精度，而采用最小步长的话，则共需1012步，这是目前的计算机技术难以负担的.当压力变化较大的时候，这一优点更加凸显，此时所需的最小步长会进一步缩小，这种情况下用固定步长法是不可能完成计算的.

4 结论

在几种用于求解Rayleigh-Plesset方程的固定步长方法中，Runge-Kutta-Fehlberg法效果最好，可以使用相对较长的步长得出准确的结果.在压力变化较小的时候，Euler法、改进Euler法也可以用于求解Rayleigh-Plesset方程，但当压力变化增大的时候，固定步长法难以得到较准确的解.为解决这个难题，提出了变步长法，即时间步长随空泡半径改变速率而自动修正的方法.

与传统的固定步长法相比，变步长法具有以下几个优点：

(1)变步长法可以模拟完整的生长溃灭循环过程，而不会在两次循环间出现奇异点;

(2)变步长法适用于较大压力变化的情况，可以得到较准确的结果;

(3)变步长法可以在保证精度的情况下，显著降低计算步数，本算例中仅为固定步长法的1/20.

因此，变步长法更适合于解Rayleigh-Plesset方程，在节约系统资源，保证结果精度等方面都有较好表现，该方法为研究Rayleigh-Plesset方程提供新的思路和算法.

［1］RAYLEIGHL.On the pressure developed in a liquid during collapse of a spherical cavity［J］.Philosophical Magazine，1917(34)：94 －98.

［2］卢义玉，葛兆龙，李晓红.空化空泡发育和溃灭过程的数值分析［J］.中国矿业大学学报，2009，38(4)：582 －594.

［3］PLESSET M S，CHAPMAN R B.Collapse of an initially spherical vapour cavity in the neighbourhood of a solid boundary［J］.Journal of Fluid Mech，1971(47)：83 －290.

［4］SHIMA A.Mechanism of impact pressure generation from spark generated bubble collapse near a wall［J］.AIAA J，1983(21)：55 －59.

［5］PLESSET M S.The dynamics of cavitation bubbles［J］.ASME Journal of Applied Mechanics，1949(16)：228 －231.

［6］王鹏.随机常微分方程数值分析中的若干方法［D］.长春：吉林大学，2008.

［7］郭晓梅.常微分方程的数值解法［J］.网络财富，2009，(4)：132.

［8］YAN Yu-bin.Smoothing a properties and approximation of time derivatives for parabolic equations：variable time steps ［J］.BIT Numerical Mathematics，2003，43(1)：647-669.

［9］PLESSET M S，PROSPERETTI A.Bubble dynamics and cavitation［J］.Annual Review Fluid Mech，1977(9)：145 －185.

［10］FUJIKAWA S，AKAMATSU T.Effects of the nonequiliblium condensation of vapour on the pressure wave produces by the collapse of a bubble in a liquid［J］.Journal of Fluid Mech，1980，97(6)：481 －512.

［11］QIN Z，BREMHORST K，ALEHOSSEIN H.Simulation of cavitation bubbles in a convergent－divergent nozzle water jet［J］.J.Fluid Mech，2007，573(4)：1－25.