雙泡耦合聲空化動力學(xué)過程模擬*

許龍 汪堯

1) (中國計量大學(xué)理學(xué)院,杭州 310018)

2) (中國計量大學(xué),浙江省智能制造質(zhì)量大數(shù)據(jù)溯源與應(yīng)用重點實驗室,杭州 310018)

為了對雙泡耦合的聲空化過程進行模擬,本文從流體動力學(xué)控制方程和流體體積分數(shù)模型出發(fā),在Fluent 軟件中構(gòu)建雙泡耦合超聲空化三維有限元仿真模型,對超聲波驅(qū)動下流體中雙泡耦合聲空化動力學(xué)過程進行數(shù)值模擬,并通過對空化氣泡周圍聲場的變化進行分析研究雙泡耦合聲空化的非線性動力學(xué)特性.結(jié)果顯示: 在超聲波驅(qū)動下,球形氣泡先緩慢擴張,擴張到最大半徑后迅速收縮直至潰滅;耦合雙氣泡間存在相互作用力,使得空化氣泡的擴張受到抑制、氣泡收縮時間增長;空化氣泡在收縮階段的能量轉(zhuǎn)換能力增強,相比單氣泡聲空化,耦合雙氣泡潰滅時氣泡內(nèi)部的壓強更大.本文分析結(jié)果將為超聲空化泡群的動力學(xué)過程模擬提供參考.

1 引言

超聲空化通常是指液體中存在的微氣泡在超聲波的作用下振蕩、生長、收縮和崩潰等一系列非線性動力學(xué)過程[1].1917 年,Rayleigh 提出Rayleigh氣泡動力學(xué)模型之后,Plesset 等在此基礎(chǔ)上提出了Rayleigh-Plesset 方程.此后大部分單泡和泡群徑向脈動的理論研究都是考慮各種復(fù)雜情況時對Rayleigh-Plesset 方程進行修改[2].文獻[3]利用高速攝影系統(tǒng)對氣泡運動和泡群的分布進行了觀察.文獻[4]對超聲波驅(qū)動下的泡群系統(tǒng)共振特性進行了研究.文獻[5]提出了在考慮氣泡耦合振動影響時球狀泡群內(nèi)氣泡振動的動力學(xué)方程.文獻[6] 對空化氣泡振動特性的影響因素進行了研究.文獻[7]研究了泡群中氣泡之間的邊界振動.文獻[8]研究了空化氣泡潰滅時間的影響因素.

以上關(guān)于空化泡的理論研究,大多數(shù)都是從空化泡進行理想的球形振蕩出發(fā),但實際研究中發(fā)現(xiàn)空化泡運動并不總是保持球形[9].2005 年,An 等[10]考慮非球形對稱聲壓擾動對空化氣泡的影響,對氣泡的不穩(wěn)定性相圖進行了計算.Liang 等[11?13]于2012 年研究了球形聲場環(huán)境中雙氣泡的非球形變形,2017 年和2021 年對非球形聲場中雙氣泡的運動特性以及基于勢流和微擾理論對超聲場中雙氣泡的不穩(wěn)定性進行了研究.馬艷等[14]對雙空化氣泡之間相互作用力對氣泡形狀的影響進行了研究.Wu 等[15]使用高速攝像機觀測了油水交界面處空化氣泡的演化過程.但是僅憑借現(xiàn)有實驗設(shè)備,很難對超聲空化過程聲場中各種變化進行精確觀測.近年來,由于計算流體動力學(xué)(computational fluid dynamics,CFD)技術(shù)的快速發(fā)展,使得可以對超聲空化過程中的空化氣泡運動過程以及聲場中各位置基本物理量(如速度、壓力、溫度等)的變化進行監(jiān)測分析[16,17].文獻[18]基于Fluent 軟件環(huán)境,模擬了近壁面空化氣泡潰滅過程.文獻[19]提出了一種運用流體體積(VOF)方法追蹤多個空化泡的運動界面的計算方法.文獻[20]采用VOF 方法模擬了管道內(nèi)空化泡的動力學(xué)特性.文獻[21]對空化效果與超聲頻率之間的關(guān)系進行了分析.文獻[22]將有限體積法和流體體積模型結(jié)合起來對近壁面微氣泡的潰滅特性進行研究.之后還有許多學(xué)者對近壁面以及靜態(tài)流場中空化泡的形變開展了研究[22].但是,對于在聲場作用下微小空化泡的形狀變化的研究卻非常少.

本文建立了雙氣泡耦合超聲空化的三維有限元仿真模型,基于Fluent 軟件對超聲波驅(qū)動下液體環(huán)境中雙氣泡耦合振動的非線性動力學(xué)過程進行模擬,研究超聲空化過程中空化氣泡形態(tài)變化與周圍物理場變化的關(guān)系.該仿真模型將為后續(xù)超聲空化泡群動力學(xué)過程模擬提供參考.

2 基本理論與仿真模型的建立

2.1 控制方程

對于超聲空化過程中的數(shù)值模擬,有如下假設(shè): 1)流體是非可壓縮牛頓流體,并且以湍流形式流動;2)忽略重力因素;3)忽略體積力作用.

根據(jù)質(zhì)量守恒定律,質(zhì)量的連續(xù)方程為[17]

式中ρ代表流體密度,μ代表流體速度.

由動量守恒定律,Navier-Stokes(N-S)方程為[17]

式中ε代表流體運動黏度,P代表流體壓力.

由于超聲空化過程中空化氣泡狀態(tài)變化涉及氣-液兩相流.所以采用氣體體積分數(shù)函數(shù)αq來追蹤流體中空化氣泡的邊界變化情況.用VOF 模型模擬氣-液交界面[23],在VOF 模型中各部分流體不能相互穿插,體積分數(shù)的連續(xù)方程及約束條件分別由(3)式和(4)式表示:

其中,q代表流體的相,αq代表流體在網(wǎng)格中所占空間的比重.

2.2 有限元模型

理論分析大多是基于球形空化氣泡的假設(shè),但實際的空化氣泡并不總是保持球形[24].為了讓數(shù)值模擬的過程更接近實際情況,構(gòu)建了圖1 所示的三維仿真分析模型,圖中的流體環(huán)境為底面半徑為0.2 mm,高0.4 mm 的圓柱體區(qū)域,初始時刻存在兩個半徑R0=0.00309 mm 的球形空化氣泡.

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

對Fluent 計算過程中的邊界狀態(tài)進行設(shè)置:1)左側(cè)邊界Pressure inlet 為聲壓入口邊界.輸入驅(qū)動聲壓函數(shù)P=Pasin(2πft) 用UDF 方式寫入;2)邊界壁面采用剛性壁面;3)空化氣泡內(nèi)部邊界設(shè)為interior.Fluent 計算過程中的參數(shù)設(shè)置如表1 所示.

表1 Fluent 參數(shù)設(shè)置Table 1.Parameter setting of Fluent.

2.3 雙泡耦合在一個聲波周期內(nèi)的超聲空化演化過程

利用Fluent 求解計算了在超聲波驅(qū)動下雙泡耦合超聲空化在一個聲波周期內(nèi)的演化過程.在整個求解過程中,將超聲波頻率設(shè)為f=25 kHz,驅(qū)動超聲波聲壓幅值設(shè)定為Pa=1.29 atm (1 atm=1.013×105Pa),由于超聲波驅(qū)動信號的周期T=40 μs,因此設(shè)置每一時間步長是0.01 μs,總步數(shù)為4000 步.圖2 所示為雙泡耦合超聲空化在一個超聲波周期內(nèi)的演化過程.

從圖2 可以看出,初始為理想球形的空化氣泡在超聲波開始輸入后先緩慢擴張,其中圖2(a)—(h)為空化氣泡擴張階段,約在t=18.80 μs 空化氣泡擴張到最大半徑;空化氣泡在擴張到最大半徑后開始收縮,約在t=22.55 μs 時潰滅.從數(shù)值模擬獲得的空化氣泡形態(tài)變化圖可以看出,在超聲波驅(qū)動下的空化氣泡形狀變化并不是規(guī)則的球形[25],而是沿聲波激勵方向逐漸往空化氣泡的兩端分裂,這也與文獻[26]實驗觀測結(jié)論相符.由于能在聲場中懸浮的空化泡都很小,一般很難直接觀測到氣泡的形狀,許多研究人員對氣泡周圍液體環(huán)境中各物理量的變化進行研究來分析氣泡的運動狀況.2001 年,Wang 等[27]為了研究非球形聲場下的氣泡動力學(xué)行為,定義了描述氣泡形變的參量,用來表征氣泡偏離球形的程度,這一參量的取值與聲場參數(shù)有關(guān).其研究表明,聲場的非球形對稱性越明顯,則氣泡的非球形對稱性越明顯.

圖2 不同時刻氣泡的形態(tài)變化情況 (a) 0 μs;(b) 4 μs;(c) 8 μs;(d) 10 μs;(e) 12 μs;(f) 14 μs;(g) 16 μs;(h) 18.80 μs;(i) 20 μs;(j) 22 μs;(k) 22.2 μs;(l) 22.4 μsFig.2.The shape change of the bubble at different time: (a) 0 μs;(b) 4 μs;(c) 8 μs;(d) 10 μs;(e) 12 μs;(f) 14 μs;(g) 16 μs;(h) 18.80 μs;(i) 20 μs;(j) 22 μs;(k) 22.2 μs;(l) 22.4 μs.

在觀察空化氣泡的演化過程中發(fā)現(xiàn),隨著空化泡體積的增大,垂直聲波激勵方向上空化氣泡兩側(cè)存在變形不對稱的現(xiàn)象,在收縮階段,這種現(xiàn)象更加明顯,為了進一步研究這一現(xiàn)象,本文對超聲空化過程中空化氣泡周圍聲場變化進行了研究.

3 雙空化泡之間的相互作用

3.1 監(jiān)測點設(shè)置

在進行超聲空化模擬過程中發(fā)現(xiàn),空化氣泡內(nèi)和氣泡周圍物理場會發(fā)生變化,于是本文在氣泡內(nèi)和氣泡外圍設(shè)置了一系列監(jiān)測點,如圖3 所示,并對超聲空化過程中這些位置壓強變化進行了數(shù)值模擬.在圖3 中,p10和p20分別表示兩個空化氣泡中心位置壓強,分別表示聲壓激勵方向上間距為10 μm 的8 個點的壓強,分別表示垂直聲壓激勵方向上間距為10 μm 的12 個點的壓強.

圖3 監(jiān)測點位置示意圖Fig.3.Schematic diagram of the monitoring points.

3.2 雙泡耦合超聲空化過程中聲場變化

本文對空化氣泡內(nèi)和氣泡周圍流場中不同位置壓強變化進行了研究,數(shù)值模擬過程中各項參數(shù)設(shè)置與表1 相同,得到了各位置的壓強變化如圖4所示.

從圖4 可以看出,在聲波激勵方向上,各監(jiān)測點壓強變化趨勢都是從標(biāo)準(zhǔn)大氣壓開始減小,降至負壓狀態(tài)后經(jīng)過一段時間開始增大至峰值,達到峰值之后在振蕩中下降至標(biāo)準(zhǔn)大氣壓;垂直于聲波激勵方向上的監(jiān)測點壓強變化趨勢與聲波激勵方向上相同;在垂直聲波激勵方向上,對比氣泡1 兩側(cè)對稱位置的壓強變化發(fā)現(xiàn),距離空化氣泡中心距離相同的位置存在空化氣泡一側(cè)的壓強數(shù)值低于無空化氣泡一側(cè),即耦合氣泡之間存在相互作用,使得雙泡之間區(qū)域的壓強減小[28].

圖4 監(jiān)測點壓強變化 (a) 聲波激勵方向;(b) 垂直聲波激勵方向Fig.4.Pressure change of monitoring points: (a) Acoustic excitation direction;(b) vertical acoustic excitation direction.

為了研究耦合空化氣泡之間的相互作用,對垂直聲波激勵方向各位置壓強進一步研究,圖5 所示是垂直聲波激勵方向上各位置壓強隨時間變化曲線.

通過觀察圖5 可以發(fā)現(xiàn),當(dāng)超聲波開始激勵時,各監(jiān)測點的壓強變化都是從標(biāo)準(zhǔn)大氣壓開始慢慢下降,這一階段空化氣泡處于緩慢擴張過程中,在這一過程中空化氣泡兩側(cè)位置對稱的監(jiān)測點(pa與ph、pb與pg、pc與pf、pd與pe)的壓強變化曲線基本重合,即在垂直聲波激勵方向上距離空化氣泡中心相同距離處兩側(cè)壓強變化基本一致;從t=20.04 μs 開始,空化氣泡進入收縮階段,各位置壓強數(shù)值隨著空化氣泡半徑減小而增大,與初始階段曲線基本重合不同,ph先于pa開始變化,同樣,pg先于pb、pf先于pc、pe先于pd,即在垂直聲波激勵方向上空化氣泡兩側(cè)對稱位置壓強變化出現(xiàn)明顯差異,存在耦合氣泡的一側(cè)壓強先于無耦合氣泡一側(cè)開始變化.

圖5 垂直聲波激勵方向監(jiān)測點壓強變化Fig.5.Pressure changes at monitoring points in the direction of vertical acoustic excitation.

為了對這一現(xiàn)象進行進一步研究,將同樣驅(qū)動聲場中只存在單氣泡和耦合雙氣泡在收縮階段壓強變化曲線進行比對,如圖6 所示.

圖6 塌縮階段監(jiān)測點壓強變化 (a) 單氣泡;(b) 耦合氣泡Fig.6.Monitor pressure changes at the point during collapse: (a) Single bubble;(b) coupled bubble.

從圖6 可以看出,在空化氣泡收縮階段,單氣泡和耦合雙氣泡附近壓強變化都是在振蕩中增大,耦合雙氣泡附近的壓強變化振蕩幅度更大;單氣泡和耦合雙氣泡兩側(cè)對稱位置壓強開始變化的時刻不同,將單氣泡和耦合雙氣泡在垂直聲波激勵方向上兩側(cè)對稱位置壓強開始變化時間差進行比對,如表2 所示.

表2 空化泡兩側(cè)位置壓強變化起始時刻差Table 2.Initial time difference between pressure changes on both sides of cavitation bubble.

從表2 可以看出,在垂直于聲波激勵方向上,空化氣泡兩側(cè)對稱位置壓強變化時間差沿著氣泡中心向外逐漸增大.以距離空化氣泡中心40 μm處位置為例,單氣泡在該位置壓強開始增大的時刻約為21.17 μs,兩側(cè)對稱位置開始增大的時間相差約為0.06 μs;雙氣泡耦合時,存在耦合氣泡的一側(cè)(兩泡之間的位置)壓強開始增大的時刻為20.02 μs,與另外一側(cè)對稱位置壓強開始增大的時間相差約1.02 μs.對比單氣泡可以發(fā)現(xiàn),耦合雙氣泡彼此相互作用,相互作用力對空化泡的徑向振動周期產(chǎn)生影響.相比于單氣泡,耦合雙氣泡擴張階段縮短,更早進入收縮階段.采用Matlab 圖像處理算法,對各個階段空化泡的形態(tài)圖像加以處理,對圖像進行計算并轉(zhuǎn)化為空化氣泡截面面積,從而計算出各個時刻空化氣泡的歸一化半徑Rmax/R0.表3 列出了相同超聲波驅(qū)動下單氣泡和雙氣泡耦合時空化氣泡擴張階段半徑的最大值、潰滅時刻空化泡內(nèi)部壓強最大值及其出現(xiàn)的時刻.

表3 不同空化氣泡最大半徑和氣泡內(nèi)部最大壓強(Pa=1.29 at m,f=25 kHz,R10=3.09 μm,R20=3.09 μm)Table 3.The maximum of radius and inside pressure of two kin ds of bubbles (Pa=1.29 atm,f=25 kHz,R10=3.09 μm,R20=3.09 μm).

從表3 可以看出: 在空化氣泡擴張階段,單氣泡在t=19.20 μs 時擴張到最大,最大半徑Rmax=30.75 μm,耦合雙氣泡在t=18.80 μs 時擴張到最大,最大半徑Rmax=26.64 μm.相對于單氣泡,耦合雙氣泡在擴張階段的最大半徑較小、擴張時間較短,由此表明耦合雙氣泡之間的相互作用抑制了氣泡的擴張.空化氣泡在收縮過程中,泡內(nèi)壓強隨著氣泡收縮不斷增大,當(dāng)空化氣泡收縮至最小半徑時氣泡潰滅,發(fā)生劇烈的空化現(xiàn)象,此時氣泡中心位置壓強達到最大值.由表3 可知,單氣泡在t=22.01 μs 時氣泡中心位置壓強達到最大值Pmax=77.21 MPa,耦合雙氣泡在t=22.55 μs時氣泡中心位置壓強達到最大值Pmax=81.81 MPa,相比于單氣泡,雙泡耦合時空化氣泡收縮階段時間延長,空化氣泡潰滅時氣泡內(nèi)部壓強增大.在一個驅(qū)動超聲波周期內(nèi),單泡的徑向振動周期T=22.01 μs,耦合雙氣泡的徑向振動周期T=22.55 μs,耦合雙氣泡之間的相互作用改變了氣泡的振動周期;耦合雙氣泡潰滅時氣泡內(nèi)部的壓強較單氣泡有所增大,即耦合雙氣泡的相互作用提高了聲空化氣泡的能量轉(zhuǎn)換效率,該結(jié)論也與文獻[28]中所敘述一致.

4 結(jié)論

本文使用Fluent 軟件構(gòu)建了雙氣泡耦合超聲空化的三維有限元仿真模型,并對超聲波驅(qū)動下水中雙氣泡耦合的空化動力學(xué)過程進行數(shù)值模擬.在此基礎(chǔ)上,研究了雙氣泡耦合聲空化動力學(xué)過程與空化氣泡內(nèi)和氣泡周圍壓強變化的關(guān)系,并與相同超聲波驅(qū)動下單氣泡空化過程做了對比分析.

1)在雙氣泡耦合超聲空化過程中,耦合雙氣泡內(nèi)和泡外附近流場壓強變化趨勢是從標(biāo)準(zhǔn)大氣壓先減小至負壓狀態(tài),一段時間后從負壓狀態(tài)開始增大達到峰值,最后在振蕩中減小至標(biāo)準(zhǔn)大氣壓;空化氣泡內(nèi)外壓強變化與氣泡體積變化成反比.

2)耦合雙氣泡之間存在相互作用力,空化氣泡之間的相互作用對空化氣泡的擴張有抑制作用,在耦合雙氣泡擴張過程中,相比于單泡,耦合雙氣泡的最大半徑減小、膨脹時間縮短;對空化氣泡的收縮有加劇作用,在氣泡收縮階段,相比于單泡,耦合氣泡收縮開始時間提前、持續(xù)時間延長;空化氣泡的徑向振動周期發(fā)生了變化.

3)耦合雙氣泡相互作用于彼此,增加了空化氣泡的能量轉(zhuǎn)換能力.在空化氣泡潰滅時,耦合氣泡內(nèi)壓強大于單氣泡內(nèi)部.空化氣泡耦合使得空化現(xiàn)象更為劇烈.

本文研究結(jié)果對研究超聲空化泡群動力學(xué)行為具有一定的參考價值,后續(xù)還將繼續(xù)進行相關(guān)研究.