壓力波作用下二維橢圓形氣泡界面演化規律研究

王宇飛,劉元清,王靜竹,*,王一偉,3,黃晨光,3

(1.中國科學院 力學研究所 流固耦合系統力學重點實驗室,北京 100190;2.中國科學院大學 工程科學學院,北京 100049;3.中國科學院大學 未來技術學院,北京 100049;4.北京宇航系統工程研究所,北京 100076)

0 引 言

空化是水動力學中重要的現象,通常發生在高速航行的回轉體、水翼、以及高速運轉的螺旋槳葉片的表面低壓區。空化相變生成大量氣泡,構成空泡群。空泡群的非穩態演化特別是潰滅現象,不僅嚴重影響航行體水動力性能,產生噪聲與振動,甚至能夠剝蝕結構表面導致破壞。水中壓力波與氣泡相互作用是空泡群潰滅的基本形式,而壓力波作用下氣泡的運動響應研究是解決空泡群潰滅問題的核心基礎。

針對壓力波與氣泡相互作用,Ding 和Gracewski[1]通過數值模擬研究了在較弱壓力波作用下氣泡的響應過程,由于表面張力相對較大,直徑100 mm 的氣泡可以保持球狀進行體積振蕩。Abe等[2]采用顯微鏡擴大觀測也發現直徑小于50 mm的氣泡,在壓力波作用下更容易保持球狀振蕩。但是,當氣泡尺寸相對較大時,在壓力波作用下,氣泡會發生非球形演化,其中,生成射流是非球形演化的最顯著特征。Dear等[3-4]通過實驗觀測的手段,分析了壓力波加載二維柱狀圓形氣泡動態響應,發現射流的方向與壓力波的傳播方向相同。Bourne和Field[5-6]分析了壓力波加載下不同形狀界面的演化特征,指出射流的方向與界面的曲率有關。隨后,國內外學者利用實驗觀測和數值模擬手段分析了壓力波與氣泡耦合問題。在數值模擬研究中,為了提高壓力波和氣/液兩相界面的捕捉,建立了許多創新的數值計算方法,如Shyue的半保守模型[7-8]、高分辨率的波前追蹤方法[9]、基于自適應特征匹配[10]、高精度算法[11-12]、自由-拉格朗日方法[13]等,分析了壓力波在界面的波系變化、界面渦量分布和射流生成等關鍵過程。但是,射流生成的內在機理尚不十分明確。

從現象上看,射流的生成與RM 不穩定現象類似,當壓力波加載不同密度的物質界面時,由于界面處壓力梯度與密度梯度不共線,導致斜壓機制,進而界面獲得加速度,界面上的擾動會逐漸演化發展[14-16]。RM 不穩定性是一類十分復雜的多尺度/強非線性問題,是高能量密度物理研究的主要內容之一,廣泛應用于慣性約束聚變、尖端武器和天體物理等領域。目前在RM 不穩定性的實驗研究中,主要利用沖擊波加載不同密度氣體的界面。本研究中的液/氣界面作為不同密度流體介質的一個特例,界面的演化現象與氣/氣界面的RM 不穩定性相似,主要體現在,當壓力波加載界面時,兩者都涉及各種流體動力學現象的強烈耦合問題,包括壓力波在界面處的反射、透射、繞射等波系變化,渦量的產生與分布。但是由于界面物理特性方面的差異,兩者也有不同。因此,水中壓力波作用氣泡界面的演化在學術研究中有著重要的研究價值與應用背景。

本研究以水中壓力波加載不同傾斜角度的二維橢圓形氣泡為主要對象,通過數值模擬方法,分析壓力波作用下氣泡界面的非球形演化過程,獲得了氣泡傾斜角度對射流生成與發展的影響規律,探索了水中壓力波作用下射流生成的內在力學機理。在數值模擬中,基于OpenFOAM 開源程序的二次開發,利用可壓縮的VOF方法和LES方法分別捕捉氣液運動界面和流場細節信息。

1 數值模擬方法

1.1 控制方程

對于水中壓力波作用下氣泡的動態響應問題的研究,本文基于Open FOAM 開源軟件,采用可壓縮的VOF(Volume of Fluid)方法捕捉氣液運動界面的演化,質量方程有如下形式:

其中αl和αg分別為水和空氣的體積分數。

將式(1)和(2)相加可得:

由于式(3)中的對流項中含有非零項?·u,故無法直接進行求解。所以將方程(1)展開后代入式(3)可得:

為了使氣液界面更加尖銳,利用了Weller提出的界面壓縮算法,引入人工對流項抵消數值耗散帶來的誤差。因此可壓縮VOF方法的質量方程如下:

其中uc為界面壓縮速度

另外,可壓縮VOF 方法的動量方程形式和混合物方程可寫為:

其中,u 混合相速度,ρ 混合相密度,τ 黏性 力,Fσ根據CSF模型給出,即Fσ＝σκ ?αl,σ 表面張力系數,κ界面曲率,h 焓,K 動能,q 熱通量。

1.2 湍流模型

為了能夠獲得兩相界面流場的細節,本文采用LES(Large Eddy Simulation)方法。基本原理是,先在選定區域內采用濾波操作,將渦分為大尺度渦和小尺度渦,對于大渦直接通過求解N-S方程得到,對于小尺度的渦通過引入亞格子模型,體現其對于大尺度渦的影響。在大渦模擬中,經過濾波操作以后的物理量以橫杠標記:

其中,Δx1、Δx2和Δx3代表網格各個方向的尺寸。對于可壓縮流動,為了避免濾波后產生的非線性應力項,使方程封閉,本文還需采用Favre濾波,操作后的物理量用波浪線表示。濾波后控制方程中產生的亞格子應力張量τSGS通過Boussinesq假設封閉:

其中εSGS為亞格子黏性耗散率。

1.3 計算域設置

圖1為本文計算域的示意圖,建立了120 mm×180 mm×2 mm 的長方體計算域,用來呈現二維幾何形狀。在長方體中心處設置了一個圓柱形氣泡。氣泡右側設定了矩形高壓區,用于模擬水中平面壓力波,計算域的上下兩個面設置為無滑移壁面。另外四個面定義為出口,水可以自由流入流出,通過采用消波邊界條件來避免壓力波的反射。坐標系的原點建立在長方體的中心,也就氣泡中心,以長度方向為x軸,寬度方向為y 軸,高度方向為z 軸。計算開始后,壓力波沿y 軸方向傳播。以壓力波陣面恰好接觸到氣泡的瞬間為0 時刻。

圖1 計算域示意圖Fig.1 Schematic of computational domain

1.4 計算方法驗證

計算域采用結構化網格劃分,為了排除網格數密度的影響,進行了網格無關性驗證。圖2展示了利用三種密度網格計算氣泡潰滅過程的對比。高、中等及低密度網格數目分別為533萬、420萬和310萬。從結果看,對于氣泡界面的演化,不同密度網格的計算結果基本一致。下文數值模擬采用的高密度網格,網格數為533萬。

圖2 不同密度網格計算結果比較。Fig.2 Comparison of bubble interface for different-type meshing

為了驗證本文計算方法的正確性與可行性,將數值模擬的結果與文獻中的實驗結果[17]進行了對比。參考Bourne和Field研究,氣泡初始直徑是12 mm,在他們的研究中,平面壓力波由于炸藥沖擊青銅板生成的,是一個矩形波,壓力峰值是260 MPa。在本文數值模擬中,利用矩形的高壓區來模擬水中平面壓力波,再根據張阿漫等研究[18],經過大量數值計算后選取與實驗結果最接近的矩形高壓區的寬度,2 mm。圖3展示了氣泡界面演化的實驗觀測與數值模擬的對比。每個圖片的時間間隔是10μs。壓力波的傳播方向是從下向上,加載氣泡后,其下表面逐漸變平,然后向氣泡中心凸起,形成高速射流,最終可以打穿氣泡壁。總體來說,氣泡的形態變化基本與實驗觀測相吻合,從而證明了數值模擬方法的可靠性。

圖3 針對水中壓力波誘導氣泡界面變形的數值模擬和實驗結果[17]對比Fig.3 Comparison between numrical simulatin and experimental observation[17]for interfacial evolution of bubble induced by underwater pressure wave

2 結果與討論

本文主要分析水中壓力波加載不同傾斜角度的二維橢圓形氣泡的界面演化規律。橢圓氣泡的傾角定義為氣泡長軸和壓力波傳播方向的夾角。在計算中,橢球形氣泡的長軸半徑設置為4 mm,短軸半徑為2 mm,模擬水中平面壓力波的矩形高壓區的寬度為2 mm,初始溫度設置為320 K。

首先研究了橢圓形氣泡傾斜角度為0°的情況。圖4展示了壓力波加載后氣泡界面形態變化。水中壓力波傳播方向從右向左。與圓形氣泡的情況(圖3所示)類似,在壓力波傳播后,與其先接觸的氣泡界面開始變形,從凹面變成平面,再逐漸向氣泡中心運動,形成射流(主射流),與壓力波傳播方向相同。隨著氣泡界面整體運動,在t＝50μs時,在氣泡長軸的另一端也形成一個射流(次射流),與射流1方向相反。在t＝55μs時,主射流擊穿氣泡壁,造成氣泡分裂。分裂后的兩部分在主射流和次射流的共同作用下,界面不斷收縮,形成類似旋渦的結構(t＝70μs)。從圖中可以看出,在水中壓力波加載橢圓形氣泡后,在氣泡長軸兩端各生成一個射流,與壓力波傳播方向相同的射流占主導作用,因此,射流角度為0°。其中射流角度為主射流運動方向與水中壓力波傳播方向的夾角。

圖5展示了傾斜角度為35°情況下氣泡動態響應。在壓力波傳播后20μs時,發現氣泡界面明顯變化的地方是橢圓形的長軸的一端,并不是與壓力最先接觸的地方。因此,發現氣泡界面的運動除了和壓力波的壓強有關,還依賴于界面的曲率。在曲率大的地方,界面運動速度也快。由此可得,射流是在氣泡長軸一端生成的,但是并不是指向氣泡中心的。在壓力波作用下,射流運動的方向與壓力波傳播方向的夾角為31.4°,因此,該工況下射流角度為31.4°。同時,在氣泡長軸的另一端又生成了一個射流(次射流)。當兩個射流相聚并撞擊后,氣泡分裂成兩部分,界面繼續變形,之后的運動變得沒有規律。

圖4 壓力波作用橢圓形氣泡的界面演化:傾斜角度0°Fig.4 Interfacial evolution of elliptical bubble induced by pressure wave:inclination angle 0°

圖5 壓力波作用橢圓形氣泡的界面演化:傾斜角度35°Fig.5 Interfacial evolution of elliptical bubble induced by pressure wave:inclination angle 35°

圖6 展示了傾斜角度為90°時氣泡界面響應情況。在水中壓力波傳播過后,隨著界面的運動,先接觸壓力波的氣泡界面不再保持平滑,形成了一些小的擾動,而是在氣泡長軸兩端氣泡界面出現了明顯的向內收縮(t＝30μs),形成兩個射流。在該工況下,主射流和次射流的速度相同,對氣泡界面變形的作用相同。他們的速度方向與壓力波傳播方向的夾角均為53.9°,因此,該工況下射流角度為53.9°。同時,運動前期生成的小擾動逐漸增大,形成了有兩個波峰的擾動。但是,與射流相比,小擾動的幅值和速度都是非常小的。隨著兩個射流繼續向氣泡內運動,在t＝55μs時,上下的射流匯聚并碰撞,造成氣泡分裂,兩個波峰變成了一個。隨后,主體氣泡界面繼續運動,形成了一個與壓力波傳播方向相反的水平射流。該射流最終穿透氣泡壁。

圖6 壓力波作用橢圓形氣泡的界面演化:傾斜角度90°Fig.6 Interfacial evolution of elliptical bubble induced by pressure wave:inclination angle 90°

為了分析壓力波作用下氣泡射流形成的內在機制,提取了不同時刻下的界面的渦量場信息,如圖7所示。圖中的黑色實線表示氣泡的輪廓,這里定義液相體積分數0.9的等值面為氣泡邊界。本文主要針對是二維幾何,所以與x 和y 軸相比,z 軸方向距離很小,所以氣泡變形和渦量主要存在x-y 平面上。圖7展示的是x-y 平面上的渦量場。在t＝25μs時刻,氣泡上部的射流是在一對渦量的作用下生成的,右邊是逆時針方向的旋渦,左邊是順時針的。同時,界面處生成的小擾動也是在渦對作用下生成的。在t＝55μs時刻之后,兩個射流相聚并穿透氣泡壁后,因此,作用在兩個射流上的兩對渦量變成了一對(如t＝55μs),上面旋渦的方向為逆時針,下面的為瞬時針。在這樣的渦對作用下,氣泡生成了與壓力波傳播方向相反的水平射流。隨著射流在泡內運動,渦量強度逐漸減弱。

為了明確誘導射流和界面擾動的渦量生成機制,提取了氣泡界面處的斜壓場(如圖8所示)。圖中的黑色實線表示氣泡的輪廓,這里定義液相體積分數0.9的等值面為氣泡邊界。從圖中可以看出,斜壓分布與渦量分布相似,在壓力波傳播后,氣泡界面處密度梯度指向形心,壓強梯度與界面曲率垂直,所以兩者在界面處不共線,進而誘導了斜壓機制,產生旋渦。從而也揭示了氣泡界面曲率大的地方,界面運動速度快。因為曲率越大的地方,壓強梯度與密度梯度夾角越大,進而產生渦量的強度越大。

圖7 壓力波誘導氣泡界面處渦量場Fig.7 Vorticity field around the surface of the bubble induced by pressure wave

圖8 壓力波誘導氣泡界面處斜壓場Fig.8 Baroclinity field around the surface of the bubble induced by pressure wave

為了進一步揭示斜壓機制在渦量中的主導作用,分析了渦量方程中其他分項的作用。渦量方程式如下所示:

其中,式中ω、u、ν、▽ρ 和▽p 分別表示渦量矢量、速度矢量、黏度系數、密度梯度和壓力梯度。式子右邊第一項就是斜壓項,(▽ρ×▽p)/ρ2,表示由于密度梯度和壓力梯度不共線時誘導的渦度。第二項渦拉伸項,(ω·▽)u,發生在速度場不均勻的情況,特別是在三維情況下發生在兩相界面。第三項渦膨脹項,ω(▽·u),主要是由于高壓沖擊強度而產生的壓縮性。第四項黏性項,ν▽2ω,表示由于流動中的黏性效應而產生的渦度。為了定量對比上述作用項對渦量的影響,將式(13)在x、y 和z 方向上展開,可以得到:

在這里,主要討論x-y 平面上的氣泡變形和渦量。

圖9展示了在氣泡響應過程中渦量方程式中各項云圖的對比。時間分別為t＝20μs、40μs、50μs和65μs。在圖中,渦量方程式中各項顏色條的范圍是不同的,前三項,斜壓項、渦拉伸項和渦膨脹項的范圍是－5×107～5×107,而黏性項是－5×106～5×106。由此可以看出,在氣泡響應過程中,黏性對渦量演化的影響最小。從圖中前三項的云圖中,可以看出斜壓項的作用最大,第二項渦拉伸項和第三項渦膨脹項的差不多,但是根據渦量方程式,第二項和第三項方向相反。由此可見,斜壓項是占主導作用的,尤其是在氣泡響應后期階段(t＝65μs)。因此,可以得出氣泡射流的生成是由于RM 不穩定性中斜壓機制導致的,壓力波傳播后,界面處的壓力場梯度與密度梯度不共線引起的。

圖9 針對氣泡響應過程中渦量方程中各項云圖對比Fig.9 Cloud chat of each term in the vorticity equation for bubble response

圖10 壓力波作用橢圓形氣泡的界面演化:傾斜角度145°Fig.10 Interfacial evolution of elliptical bubble induced by pressure wave:inclination angle 145°

圖10 展示了傾斜角度等于145°情況下,氣泡界面演化。相比傾斜角度35°工況(圖5所示)結果,兩者的氣泡界面演化在垂直面上對稱。在壓力波掃過氣泡后,在氣泡長軸上端點位置開始出現明顯向內收縮,生成主射流,其運動方向與壓力波傳播方向夾角為31.4°,因此,該工況下射流角度為31.4°。

最后,分析了不同氣泡傾斜角度對射流角度的影響,如圖所示11。橫軸是氣泡傾斜角度,縱軸是射流角度。氣泡傾斜角度為長軸與壓力波傳播方向之間的夾角,射流角度為主射流的運動方向與壓力波傳播方向之間的夾角。當氣泡長軸與壓力波傳播方向平行的情況,射流運動方向與其相同。當氣泡長軸垂直于壓力波傳播方向的情況,射流的角度等于53.9°。氣泡傾斜角度在0°～90°的范圍,射流角度隨傾斜角度增大而增加。在90°～180°之間范圍,射流角度隨傾斜角度增大而減小。而且,從圖中還發現,射流角度在這兩個范圍內是對稱的。

圖11 氣泡傾斜角度與射流角度的關系Fig.11 Relation between inclination angle of elliptical bubble and jet angle

3 結 論

本文基于OpenFOAM 開源程序的二次開發,利用可壓縮的VOF方法和LES方法分別捕捉氣/液運動界面和流場細節信息。以水中平面壓力波加載不同傾斜角度的二維橢圓形性氣泡為主要研究對象,分析了氣泡界面演化規律,揭示了射流生成的內在機理。獲得的主要結論如下:

1)射流生成位置與氣泡傾斜角度沒有關系,都是從氣泡長軸兩個端點上生成的,然后向氣泡內運動。這可能與界面曲率有關,長軸端點的界面曲率大,因而界面運動速度快;

2)射流角度與氣泡傾斜角度密切相關。當氣泡傾斜角度等于0°的情況,主射流的運動方向與壓力波傳播方向相同,當氣傾斜角度等于90°的情況,射流的運動方向與壓力波傳播方向夾角為53.9°,當傾斜角度在0°～90°之間,射流角度隨著氣泡傾斜角度增大而增大。而且還發現,射流角度在傾斜角度90°～180°的工況與0～90°的工況是對稱的;

3)通過定量對比渦量方程中斜壓項與其他項作用,發現斜壓是誘導渦量產生的主導機制。當壓力波掃過氣泡后,氣泡界面處的壓力梯度和密度梯度不共線,誘導了斜壓機制,進而生成射流和界面的擾動。另一方面也解釋了界面曲率大的地方運動速度快,曲率越大,界面處的壓強梯度和密度梯度方向差別越大,斜壓項作用越大,從而誘導的渦量強度越大。