空化水射流空泡潰滅過程的數(shù)值分析

楊博凱，盧義玉，楊曉峰，馮明濤，張 賽

(重慶大學(xué)資源及環(huán)境科學(xué)學(xué)院，重慶400030)

0 引言

目前，高壓空化水射流已經(jīng)應(yīng)用于很多領(lǐng)域，如巖石破碎、海上石油鉆井、船舶清洗、空化水射流處理有機物廢水以及醫(yī)療技術(shù)等，但目前對于空泡生成和潰滅的機理還不清楚.

對于空化水射流而言，其空泡半徑與壓力、時間存在確定的函數(shù)關(guān)系，而三者的關(guān)系直接影響空化射流的效率以及能量釋放程度.因此，旨在研究它們之間關(guān)系的空泡動力學(xué)方程在空化理論中占有重要的地位.

對于一個包含空泡的液體，如果已知射流中單個空泡的內(nèi)外壓力，就可以計算它的半徑R［1］，同時能夠建立一個Rayleigh-Plesset方程來近似求得氣泡的運動狀態(tài)［2］.理論上這個方程可以使用常見的任意一種有限元差分法解決，但由于空泡半徑隨時間變化的過程是非線性的，常用的迭代方法在求解Rayleigh-Plesset方程時存在困難，計算量大、精確度低，特別是在空泡潰滅階段會產(chǎn)生振蕩.筆者根據(jù)變步長法的基本思想，結(jié)合自主研發(fā)的空化水射流對各種迭代法的計算效率進(jìn)行分析，提出了變步長的Rayleigh-Plesset方程計算方法.變步長法能夠更有效地模擬空泡最終潰滅階段的能量釋放過程，而這個能量釋放過程正是研究其對附近邊界產(chǎn)生破壞作用機理的關(guān)鍵［3－4］.

1 Rayleigh-Plesset空泡運動方程

Rayleigh-Plesset方程是解決空泡運動方程的經(jīng)典方程，在空化水射流的研究中，經(jīng)常需要確定空泡隨時間的變化以及最終潰滅放能情況，而通過Rayleigh-Plesset方程可以較為準(zhǔn)確地得到空泡半徑與時間的關(guān)系.它是一個二階微分方程［5］.

式中：ρ為液體密度;μ為液體運動黏度;t為時間;p∞為液體靜壓;pB為空泡內(nèi)的壓力.

式(1)可以轉(zhuǎn)換為以下的一階方程

當(dāng)x為因變量時，令R=x，可同時求解式(2)和(3)以得到方程的解.令：

采用 Euler法、Central法、改進(jìn) Euler法和Runge－Kutta－Fehlberg法4種有限元拆分法求解［6－7］，則需將方程(1)改寫成

先將式(7)轉(zhuǎn)變成形如式(2)和式(3)的兩個一階等式，以便應(yīng)用Runge－Kutta－Fehlberg法，令 x=R，則

假定 pB和 p∞以及初始條件已知，便可以通過式(8)和(9)求解空泡半徑R，即空泡從生長到潰滅的過程中的大小可以以時間的函數(shù)來表示.

在式(9)中，變量是pB和p∞.推測空泡中含有雜質(zhì)氣體，其壓力為pgo大小為Ro，水蒸氣壓力為pv，在氣體不會凝結(jié)的前提下，它是一個與溫度相關(guān)的定值［9－10］，因此，筆者將其簡化為

則空泡壓力可以表示為

液體壓力p∞(t)是一個重要的參數(shù)，它關(guān)系到空泡成長和潰滅的關(guān)鍵參數(shù)p∞.沒有可以直接求得p∞的準(zhǔn)確表達(dá)式，一般由經(jīng)驗公式來獲得，p∞由CFD軟件模擬的結(jié)果獲得.

為解Rayleigh-Plesset方程，需要先確定某些參數(shù).筆者的測試條件為產(chǎn)生空泡的介質(zhì)為水，且溫度保持在300 K不變，水的性能見表1.

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

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

2.1 固定步長法的基本算法

在滿足前述解Rayleigh-Plesset方程各項假設(shè)的前提下，對于固定步長法的計算方法可表示如下：(1)令 h0=Δt;(2)取 ti=ti－1+h0i;(3)將結(jié)果帶入式(8)和(9)中;(4)由Rayleigh-Plesset方程求得時刻的空泡半徑.

2.2 幾種數(shù)值解法的穩(wěn)定性分析

Rayleigh-Plesset方程是一個非線性微分方程，且沒有解析解，不易比較判斷各種數(shù)值方法的可靠性，因此采用下式來進(jìn)行測試.

式(12)在x=90°時具有奇異性，與原式性質(zhì)相同，因此符合要求［8］.應(yīng)用式(12)并利用 EXCEL工具，得到的4種數(shù)值方法從Δx=1°開始的步長算法分析如圖1.

圖1 對比分析幾種方程解的情況Fig.1 Comparison of solutions from different schemes with analytical solution

如圖1所示，在x不接近90°時，各種方法都能得到較為準(zhǔn)確的結(jié)果，但在整個0°≤x≤90°的范圍中考慮的話，Runge-Kutta-Fehlberg法效果最好.

4種方法中只有Runge-Kutta-Fehlberg法可以估算誤差.對比解析解可以看出Runge-Kutta-Fehlberg法的估算誤差較為保守一些.

2.3 較小壓力下各種解法解Rayleigh-Plesset方程

如表2所示，壓力從12 kPa線性下降到最低值－1 kPa，之后又恢復(fù)到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法在空泡潰滅再循環(huán)的時候產(chǎn)生嚴(yán)重的數(shù)值振蕩之外，其他三種方法得到了穩(wěn)定且?guī)缀跏窍嗤慕?

2.4 壓力變化較大時解Rayleigh-Plesset方程

根據(jù)K.BREMHORST等人的研究可知，在較大壓力變化的情況下，所需的時間步長非常小.因此，擴(kuò)大壓力的變化范圍進(jìn)行實驗，使最高和最低壓力分別為120 kPa和－10 kPa，如表3所示.當(dāng)時間步長為2×10－9s時，即使在一個生長潰滅循環(huán)之內(nèi)，4種方法也都不能完成.

表3 擴(kuò)大壓力變化時的測試情況Tab.3 Enlarged pressure variation for testing numerical methods

為盡量完成至少一個循環(huán)，縮減時間步長至1×20－9s再次計算，但從結(jié)果上看不出有明顯的變化.這是因為當(dāng)空泡即將潰滅時，時間步長與空泡半徑達(dá)到最小值的變化率相比過大，所以除非大幅縮減時間步長，否則不會有明顯的改觀.

4種方法在潰滅結(jié)束、開始第二次循環(huán)的時候都出現(xiàn)了振蕩現(xiàn)象，沒有給出合適的解.這表明，由于潰滅的最后階段壓力變化的速度和空泡半徑的變化速度都很快，固定時間步長的方法解Rayleigh-Plesset方程需要一個極小的時間步長，而過小的步長將導(dǎo)致計算步數(shù)與計算時間成幾何倍數(shù)增長，這種增長往往超過了現(xiàn)有計算機所能承受的極限.因此，必須找到一種更為合理的算法，才能在壓力變化比較大的情況下較好的解決Rayleigh-Plesset方程.

3 變步長法解Rayleigh-Plesset方程

3.1 變步長法的基本算法

在滿足前述解Rayleigh-Plesset方程各項假設(shè)的前提下，對于變步長法的計算方法可表示如下：(1)令 h0= Δt，hi=hi－1Ri/Ri－1;(2)取 ti=ti－1+hi;(3)將結(jié)果帶入式(8)和式(9)中;(4)由Rayleigh-Plesset方程求得ti時刻的空泡半徑Ri.

3.2 可變步長法及其在壓力變化較大時的測試

Rayleigh-Plesset方程的特點是可以描述空泡生長至潰滅的全過程，當(dāng)空泡半徑很小時，它的導(dǎo)數(shù)變得很大，因此需要非常小的時間步長.從上述計算情況來看，對于二階Rayleigh-Plesset微分方程，選取的時間步長并不合適.

基于上述考慮，提出了一種時間步長隨空泡半徑變化速率Ri/Ri－1而自動調(diào)整的算法，配合一些簡單的算法，如Euler法，來達(dá)到變步長解Rayleigh-Plesset方程的目的.圖5顯示了使用這種變步長法在壓力變化較大時測得的結(jié)果，與從圖2到圖4的結(jié)果相比，這無疑是一個更為準(zhǔn)確的結(jié)果.更重要的是，它在的時候并不會產(chǎn)生奇異點，因此可以連續(xù)地表示一個循環(huán)結(jié)束到下一循環(huán)開始的情況.使用這種變步長結(jié)合Euler法，計算空泡4次生長潰滅循環(huán)僅需不到1 500步.圖6是空泡半徑R隨時間變化的情況.圖中顯示，空泡結(jié)束第一次循環(huán)，完成潰滅并開始反彈的時間是t=0.000 015 s，這與固定步長法測得的結(jié)果相同，但固定步長法無法模擬潰滅結(jié)束并開始反彈的過程.

3.3 變步長法與固定步長法的比較

為對比變步長法與固定步長法的計算效果，使用固定步長法中較精確的Runge-Kutta-Fehlberg法與變步長法進(jìn)行對比.

在較小壓力變化的情況下，分別采用變步長配合Euler法與Runge-Kutta-Fehlberg法對一個生長潰滅循環(huán)計算，其中 pmax=12 kPa，pmin=－1 kPa.結(jié)果見圖 7.

如圖所示，兩種方法的結(jié)果沒有顯著的差異，均可模擬在壓力變化不大情況下潰滅和反彈的過程.然而在迭代步數(shù)上卻有很大差異.Runge-Kutta-Fehlberg法共使用16 000步，而變步長配合Euler法僅使用800步就得到了正確的結(jié)果.

圖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.若使用固定步長法，則必須采用最小的步長以達(dá)到相同的精度，而采用最小步長的話，則共需1012步，這是目前的計算機技術(shù)難以負(fù)擔(dān)的.當(dāng)壓力變化較大的時候，這一優(yōu)點更加凸顯，此時所需的最小步長會進(jìn)一步縮小，這種情況下用固定步長法是不可能完成計算的.

4 結(jié)論

在幾種用于求解Rayleigh-Plesset方程的固定步長方法中，Runge-Kutta-Fehlberg法效果最好，可以使用相對較長的步長得出準(zhǔn)確的結(jié)果.在壓力變化較小的時候，Euler法、改進(jìn)Euler法也可以用于求解Rayleigh-Plesset方程，但當(dāng)壓力變化增大的時候，固定步長法難以得到較準(zhǔn)確的解.為解決這個難題，提出了變步長法，即時間步長隨空泡半徑改變速率而自動修正的方法.

與傳統(tǒng)的固定步長法相比，變步長法具有以下幾個優(yōu)點：

(1)變步長法可以模擬完整的生長潰滅循環(huán)過程，而不會在兩次循環(huán)間出現(xiàn)奇異點;

(2)變步長法適用于較大壓力變化的情況，可以得到較準(zhǔn)確的結(jié)果;

(3)變步長法可以在保證精度的情況下，顯著降低計算步數(shù)，本算例中僅為固定步長法的1/20.

因此，變步長法更適合于解Rayleigh-Plesset方程，在節(jié)約系統(tǒng)資源，保證結(jié)果精度等方面都有較好表現(xiàn)，該方法為研究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］盧義玉，葛兆龍，李曉紅.空化空泡發(fā)育和潰滅過程的數(shù)值分析［J］.中國礦業(yè)大學(xué)學(xué)報，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］王鵬.隨機常微分方程數(shù)值分析中的若干方法［D］.長春：吉林大學(xué)，2008.

［7］郭曉梅.常微分方程的數(shù)值解法［J］.網(wǎng)絡(luò)財富，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.