單泡超聲空化仿真模型的建立及其動力學過程模擬*

徐 珂 許 龍

(中國計量大學理學院 杭州 310018)

0 引言

超聲空化通常是指液體中存在的微氣泡在超聲波作用下振蕩、生長、收縮和崩潰等一系列非線性動力學過程[1]。1895年，Thorneycroft 等[2]在觀察到經過短時間航行的潛艇螺旋槳推進器的空蝕現(xiàn)象之后，第一次提出了空化的概念。1917年，Rayleigh[3]給出了Rayleigh 氣泡動力學模型，Plesset[4]、Noltingk 等[5]、Neppiras 等[6]和Poritsky[7]考慮了液體的表面張力、黏滯性和驅動聲場等因素，給出Rayleigh-Plesset 方程。后面關于泡及泡群的理論研究[8-14]主要基于球形空化泡在考慮更復雜情況下的Rayleigh-Plesset方程的修正。

以上理論研究都是在基于空化泡球形變化基礎上做的分析，實際情況下的空化泡運動過程并不會一直保持球形。為了觀察實際情況下的動力學過程，許多學者進行了大量實驗研究[15-22]。但實驗研究受實驗設備的限制很難準確觀察流體中各位置相關物理參數(shù)。隨著計算流體動力學(Computational fluid dynamics，CFD)技術的發(fā)展，觀察空化泡運動過程中空化泡形態(tài)以及流體各位置基本物理量成為可能。目前基于CFD 的空化泡的仿真模擬研究，多數(shù)工作[23-29]都是對近壁面以及靜態(tài)流場中半徑比較大的空化泡的潰滅形變研究，很少有人研究聲場作用下微小空化泡的形態(tài)變化。

本文基于有限元仿真方法，構建了流體環(huán)境中單泡的有限元仿真模型，采用流體動力學控制方程和流體體積分數(shù)(Volume of fluid，VOF)模型，模擬了超聲驅動下水中單泡的演化過程并與實驗擬合結果進行對比，進一步探討了單泡形態(tài)變化過程中泡內及泡外單泡附近壓強和氣體密度分布。用有限元軟件對單泡的動力學行為進行仿真，仿真分析結果與實驗結果符合良好。本文建立的單泡超聲空化有限元仿真模型為后續(xù)雙泡以及泡群超聲空化動力學過程模擬提供有益參考。

1 仿真模型

1.1 控制方程

在計算空化泡演變過程中引入基本假設條件：(1)液體為不可壓縮牛頓流體，且為湍流流動；(2)不考慮重力作用影響；(3)忽略體積力的作用。

質量連續(xù)方程為[30]

式(1)中，ρ為流體密度，μ為流體速度。

Navier-Stokes(N-S)方程為[30]

式(2)中，ε為流體運動黏度，P表示流體壓力。

空化氣泡狀態(tài)變化整個過程涉及到氣-液兩相流。為追蹤水中氣泡的邊界變化，引入體積分率函數(shù)αq，采用VOF[31]進行氣-液交界面的模擬。VOF模型中流體沒有相互穿插，體積分數(shù)的連續(xù)方程(3)及約束條件(4)如下：

其中，q為流體的相，體積分率函數(shù)αq表示流體在網格中所占空間的比例。

1.2 有限元模型

為了仿真超聲驅動下水中單氣泡的演化過程，建立了如圖1所示的二維軸對稱模型，圖中流體區(qū)域為長0.6 mm、寬0.3 mm的矩形區(qū)域，其中半圓形區(qū)域為初始半徑為R的空泡。

圖1 計算模型Fig.1 Model for the calculation

計算區(qū)域中的5 個邊界條件分別為：(1)左邊axis為軸對稱邊界條件；(2)右邊界和下邊界wall為剛性壁面；(3)空泡界面interior為內部邊界；(4)上邊pressure inlet為聲壓入口邊界。為了模擬壓力入口超聲變化，使用UDF方法編寫輸入驅動聲壓函數(shù)P=Pasin(2πft)，其中Pa為驅動聲壓幅值，f為超聲頻率。

采用有限元軟件求解器中基于壓力的求解器進行計算，過程為瞬態(tài)，相關參數(shù)設置如表1所示。

表1 有限元法相關參數(shù)設置Table 1 Parameter setting of finite element method

1.3 仿真模型驗證

為了驗證本文建立的有限元仿真模型的可靠性，依據文獻[19]中超聲空化的相關實驗數(shù)據，對初始半徑R0= 6.18 μm 的超聲空化泡的動力學行為進行了仿真模擬。根據文獻[19]，計算過程中驅動聲壓幅值Pa= 1.29 atm (1 atm = 101325 Pa)，頻率f= 25 kHz，根據超聲頻率可知驅動信號的周期為40 μs，設置時間步長為0.02 μs，步數(shù)為2000 步。仿真計算的空化泡的形態(tài)變化過程如圖2所示。

圖2(a)～(j)為氣泡膨脹階段，約占驅動信號的二分之一周期，其中圖2(j)為單泡最大半徑時圖像；圖2(j)～(n)為氣泡塌縮階段，約占驅動信號的十分之一周期，21 μs以后為單泡膨脹壓縮做往復運動階段。由此可以看出氣泡形態(tài)隨時間的演化規(guī)律是先緩慢膨脹，到達最大值后，迅速塌縮后消失。由仿真計算的空化泡的形態(tài)變化圖可知，球形氣泡在超聲波驅動下其空化泡的變化形狀并非理想的球形，而是沿聲壓激勵方向向兩邊分裂。這也與相關實驗觀察[20-21]所得結論一致。

圖2 不同時刻氣泡的形態(tài)變化情況Fig.2 The shape change of the bubble at different time

基于上述超聲驅動下氣泡演變過程的有限元仿真模型，通過計算獲得不同時刻空化泡的形態(tài)圖像后，用圖像處理算法，通過讀取不同時刻的氣泡圖像，將氣泡區(qū)域的像素進行計數(shù)并轉換成氣泡面積，根據氣泡面積計算出不同時刻的氣泡半徑，從而獲得如圖3藍色點畫線所示的氣泡歸一化半徑隨時間的演化曲線。為了驗證本文仿真模型的可靠性，將本文仿真計算結果與文獻[19]中用實驗測試數(shù)據擬合的R-P方程的計算結果(圖3中的紅色點畫線)進行對比。

圖3 有限元仿真與實驗擬合的R-P 方程計算曲線Fig.3 The calculation curve of R-P equation fitted with the experiment and Finite element simulation

由圖3可知，仿真計算所得空化泡的膨脹與潰滅曲線與經實驗測試數(shù)據擬合的R-P 方程計算曲線基本一致。不同之處在于：(1)仿真結果中空化泡無法觀測到反彈階段，這是因為單泡初始半徑過小，在微泡縮小時無法適應網格大小；(2)實驗擬合的R-P 方程計算曲線略高于仿真曲線，其原因可能是仿真計算過程中考慮了湍流，而用R-P 方程計算時未考慮湍流影響，從而帶來二者的差別。

1.4 監(jiān)測點設置

單泡空化過程中泡內和泡外附近流場中各物理量會發(fā)生劇烈變化。在仿真計算過程中，為了監(jiān)測超聲激發(fā)單泡空化過程中泡內及泡外流體域內壓強和密度隨時間的變化，在泡內和泡外設置了9 個監(jiān)測點，如圖4所示。點p0位于單泡中心，點p1、p2、p3、p4在聲壓激勵方向，分別距離單泡中心10 μm、20 μm、30 μm 和40 μm，點p5、p6、p7、p8在聲壓垂直方向，分別距離單泡中心10 μm、20 μm、30 μm和40 μm。

圖4 監(jiān)測點位置示意圖Fig.4 Schematic diagram of monitoring point location

2 結果與討論

2.1 超聲作用下單泡泡內及附近壓強變化

對單泡泡內和單泡附近流域壓強變化進行了研究。仿真過程中單泡初始半徑R0= 6.18 μm，驅動聲壓幅值Pa= 1.29 atm、頻率f= 25 kHz。圖5為本文圖4中9個監(jiān)測點處壓強隨時間變化曲線。

圖5 監(jiān)測點壓強變化Fig.5 Pressure change of monitoring point

如圖5中p0曲線所示，單泡泡內壓強先緩慢減小后迅速增大，在t= 20.60 μs 達到最大值87 atm后迅速減小并做往復變化。結合圖2和圖3發(fā)現(xiàn)，泡內壓強變化與單泡體積變化成反比。

當t小于16.30 μs 時，單泡處于膨脹階段。如圖5所示，在t= 6.40 μs 之前，p1與p5曲線、p2與p6曲線、p3與p7曲線和p4與p8曲線分別重合且逐漸遞減，表明在單泡膨脹階段，泡外壓強在聲波激勵方向和聲波垂直方向上與泡中心相同距離處壓強相同，兩個方向的泡外壓強隨著與泡中心距離的增加逐漸遞減。當t= 6.40 μs，單泡半徑達到R= 10 μm，此時單泡的表面正好膨脹到p1和p5監(jiān)測點的位置，隨著時間的推移單泡進一步膨脹，p1、p5進入泡內，p1和p5曲線加速向p0曲線靠攏，并在7.2 μs時與p0曲線重合，從此刻開始p1、p5和p0點的壓強相同。在t= 10.6 μs 時，單泡半徑R=20 μm，此時p2和p6監(jiān)測點也正好位于單泡表面處，隨著時間的增加單泡進一步膨脹，p2和p6監(jiān)測點也進入泡內，p2、p6曲線也開始加速向p0曲線靠攏，并在12 μs 時與p0曲線重合，從此刻開始p2、p6與p1、p5和p0點的壓強都相同。由此可見，在單泡膨脹過程中，當監(jiān)測點剛進入泡內時，相關檢測點的壓強與泡中心的壓強有一定差別，但隨著時間的推移，相關監(jiān)測點處壓強會加速與泡中心點壓強趨于一致。

當16.30 μs＜t ＜20.60 μs 時，泡外壓強大于泡內壓強，單泡處于壓縮階段，隨著泡的壓縮，泡內壓強開始增加。在單泡壓縮階段，泡外聲壓在聲波激勵方向和垂直方向壓強變化出現(xiàn)明顯差異，從而導致在超聲空化過程中氣泡呈非規(guī)則球形變化。

當20.60 μs＜t ＜40 μs 時，如圖5所示，泡外壓強隨距離單泡中心點的距離增大而減小。同時，在聲壓垂直方向的壓強值要高于聲壓激勵方向。對比單泡形態(tài)變化，就是在單泡壓縮階段以后，單泡不再是規(guī)則的球形，并且沿聲壓激勵方向向兩邊分裂[20]，導致聲壓垂直方向的壓強值要高于聲壓激勵方向。

2.2 超聲作用下單泡泡內及附近氣體密度變化

在與上述相同邊界條件下對單泡泡內及泡外附近氣體密度變化進行了研究。圖6為圖4中9 個監(jiān)測位置處氣體密度隨時間變化曲線。

如圖6中p0曲線所示，單泡泡內氣體密度先緩慢減小后迅速增大，在t= 20.60 μs 達到最大值103 kg/m3，然后迅速減小并做往復變化。結合圖2和圖3發(fā)現(xiàn)，泡內氣體密度變化與單泡體積變化成反比。由圖6可知，單泡在緩慢膨脹階段，泡內密度各處相同且逐漸減小，泡外各處密度隨距泡心距離的增大而減小，且泡外密度在聲壓激勵方向和聲壓垂直方向始終相同。單泡在壓縮階段，單泡附近氣體密度增大，并且聲壓激勵方向氣體密度變化相對垂直方向較大，在聲壓垂直方向氣體密度要大于聲壓激勵方向氣體密度，對比單泡形變，說明單泡壓縮過程中發(fā)生了沿聲壓激勵方向的形變，單泡在壓縮階段沿聲壓激勵方向向兩邊分裂，與文獻[20]實驗觀測結果吻合。

圖6 監(jiān)測點的氣體密度變化曲線Fig.6 Variation curve of gas density at monitoring points

3 結論

建立了超聲空化泡的有限元仿真模型，利用有限元軟件模擬了超聲驅動下水中單泡的空化動力學過程?；诜抡婺M獲得的不同時刻空化泡的演化圖像，采用圖像處理算法，獲得了氣泡半徑隨時間的變化曲線。在此基礎上，分析了單泡形態(tài)變化過程中泡內及泡外附近壓強和氣體密度分布。結論如下：

(1)仿真計算所得氣泡半徑隨時間的演化規(guī)律是先緩慢膨脹到最大后迅速塌縮，與實驗擬合的R-P 方程計算所得氣泡半徑隨時間的變化趨勢一致，從而驗證了有限元仿真模型的可靠性。

(2)單泡形變過程中，泡內壓強與氣體密度變化與單泡體積變化成反比。

(3)單泡膨脹階段，單泡泡內壓強與氣體密度大于泡外附近壓強與氣體密度，并且壓強與氣體密度以單泡中心為球心向外遞減。

(4)單泡壓縮階段，單泡附近的壓強與氣體密度不再是球形均勻變化，而是在聲壓垂直方向要大于聲壓激勵方向，從而導致在超聲空化過程中氣泡呈非規(guī)則球形變化。

(5)單泡壓縮階段后，在聲壓垂直方向壓強與氣體密度要大于聲壓激勵方向氣體密度，對比單泡形變，說明單泡壓縮階段過后，沿聲壓激勵方向向兩邊分裂。

本文研究結果將為模擬更加復雜的超聲空化泡及泡群動力學過程提供參考借鑒，相關研究工作將在后續(xù)進行。