余 俊，劉國振，汪 俊，王海坤

(中國船舶科學研究中心，無錫 214082)

1 引 言

水下爆炸可以劃分為三個重要階段，包括裝藥的爆轟、沖擊波向外傳播和氣泡脈動等過程。水下爆炸過程中除了會向周圍流體傳遞沖擊波載荷外，爆轟產物形成的爆炸氣泡在脈動過程中也會向外傳遞沖擊載荷[1]。當爆炸氣泡附近存在結構和自由面等邊界條件時，在爆炸氣泡潰滅過程中容易產生指向結構或者背向自由面的水射流現(xiàn)象。水射流產生過程中涉及到爆炸氣泡、自由面或者結構邊界和水等多物質和多界面之間的復雜耦合過程，研究難度非常大。目前對于氣泡射流的理論分析僅限于簡單規(guī)則的情形[2-4]，而實驗研究的可重復性和可變性較差[5-7]，所以對氣泡的數(shù)值模擬研究逐漸得到重視。目前對于水下爆炸氣泡運動模擬的數(shù)值方法也比較多，其中邊界元方法(BEM)普遍應用在氣泡動力學模擬計算領域。BEM法計算維數(shù)降低，計算效率較高，在模擬三維復雜幾何形狀邊界面時具有一定的優(yōu)勢。文獻[8-10]在BEM方法的理論模型、計算方法和切割處理技術等多個方面開展了研究，使BEM方法得以不斷補充和完善。BEM法在氣泡脈動周期以及氣泡最大半徑等模擬方面具有較好的精度，但是在氣泡潰滅階段，由于水射流的高速運動以及水錘效應的存在，流場中的速度梯度較大，故需要考慮流體的可壓縮性以及大變形等影響。此時BEM法的求解精度明顯下降，難以適應高精度模擬的要求。除了邊界元方法外，氣-液兩相流模型也是模擬氣泡的重要方法，根據(jù)界面處理方法不同可以劃分為界面追蹤和界面捕捉兩大類。界面追蹤法是通過粒子點或標記等方法顯式地追蹤界面運動，比較適合簡單的界面運動情況，難以描述界面復雜運動。而界面捕捉方法則是隱式地捕捉界面位置，主要以水平集法(level-set)[11]和體積分數(shù)法(vof)[12]為代表，在復雜界面運動模擬方面取得了大量廣泛的應用。對于level-set方法，由于在數(shù)值通量守恒方面存在一定的不足，難以長時間精確捕捉界面運動過程。而vof方法是通過體積分數(shù)的對流方程來實現(xiàn)界面的捕捉，在長時間計算后界面容易出現(xiàn)彌散現(xiàn)象，界面出現(xiàn)模糊，從而影響計算精度。

由于上述計算方法在處理水下爆炸氣泡射流階段存在的困難和挑戰(zhàn)，故難以精確獲得水射流載荷的特征形式，無法對水射流載荷的沖擊損傷效應進行準確評估。為了克服上述困難，本文擬在多相可壓縮流體的 five -equation 計算模型基礎上，引入界面函數(shù)和界面密度的壓縮技術來提高水射流階段氣-液界面的捕捉精度，抑制界面彌散效應，從而建立爆炸氣泡運動的數(shù)值計算模型。該模型采用5階WENO重構與HLLC近似Riemann求解器進行空間離散，時間離散采用3階TVD Runge -Kutta 法，能夠捕捉沖擊波的傳播以及多相界面運動。針對該計算模型，首先通過三個典型的多相流問題算例進行考核，最后將該模型應用到水下近壁面爆炸過程中水射流現(xiàn)象的模擬，獲得了水射流載荷的典型特征，為水射流的產生機理及其損傷效應的深入研究提供了技術支撐。

2 計算模型

2.1 控制方程

對于水下爆炸氣泡運動等瞬態(tài)響應過程，可以忽略傳熱、粘性以及化學反應等影響，采用多相可壓縮流體的 five -equation 模型對爆炸氣泡及其周圍流體之間的相互作用過程進行描述，其二維控制方程可表示為[13]

(1)

式中Q,F和G分別為

(2)

(z1)t+u(z1)x+v(z1)y=0

(3)

流體的狀態(tài)方程采用剛性氣體狀態(tài)方程

(4)

2.2 界面壓縮與密度修正

基于體積分數(shù)的界面捕捉模型在計算過程中隨著計算時間的推進，兩相界面附近不可避免地會出現(xiàn)彌散效應，進而導致界面模糊不清。在流體產生渦旋和水錘現(xiàn)象的模擬過程中將顯著影響流場的速度、壓力和密度等關鍵物理量的計算精度，導致計算結果不可信。嚴重情況下會引起計算崩潰，從而使得計算模型的適用范圍縮小，計算能力下降。特別是針對水下爆炸氣泡脈動過程的模擬，氣泡運動周期相對沖擊波作用階段長得多，采用顯式計算格式的時間步長非常小，要想模擬氣泡潰滅產生水射流的過程所需的計算步數(shù)特別長(一般在幾萬步以上)。如果不對界面彌散效應進行有效抑制，則計算結果將嚴重失真。

圖1所示為未采用壓縮技術計算的100 kg TNT在水深400 m工況下近壁面爆炸氣泡收縮過程中崩潰形成的水射流現(xiàn)象，其中左側為固壁條件，其他為透射邊界，計算采用二維軸對稱模型。圖2為相同工況下采用壓縮技術的計算結果，圖3為上述兩種計算工況下壁面中心點水射流壓力的時程曲線。由上可知，如果不引入壓縮技術，數(shù)值計算產生的氣-液界面彌散效應將使得氣泡射流撞擊時頭部水的密度急劇降低，導致射流的水錘效應不明顯，水射流壓力降低。由圖3可以看出，水射流的峰值壓力由77.8 MPa下降到49.2 MPa，下降幅度約為36.8%。

圖2引入的界面壓縮技術借鑒了文獻[15-17]發(fā)展的界面壓縮技術。首先將單相介質的質量守恒方程表示為

(5)

式中K=εh|z1|-z1(1-z1)，μ0為引入的系數(shù)，當μ0→∞時，要保證K=0，即在單相流體單元中式(5)右端趨于0，即此時界面不是0就是1，無需進行壓縮。只有當z1∈(0,1)時，才對界面函數(shù)進行壓縮。

圖1 未帶壓縮下流場的密度和壓力云圖

圖2 帶壓縮下流場的密度和壓力云圖

在兩相流中，除了界面位置會產生彌散外，密度在界面處也會隨著時間彌散。本文直接給出質量守恒方程的界面壓縮方程為

(6)

式中H(z1)=tanh{[100z1(1-z1)]2}。

圖3 壁面中心點處水射流壓力曲線對比

3 數(shù)值離散

由于計算域中載荷強間斷以及接觸間斷面的存在，采用有限體積格式進行空間離散[18,19]。將控制方程(1)在控制單元Ii j上進行平均積分，化簡為

(7)

4 數(shù)值測試與驗證

4.1 一維激波管問題

為了考核計算模型對沖擊波傳播以及界面運動的捕捉能力，采用一維Sod問題來進行測試，模擬的是較為復雜的Sod激波管問題，一般是用來分析水下爆炸的氣-液兩相運動問題[22,23]，無量綱參數(shù)的初始條件為

(8)

4.2 二維激波管問題

圖4 計算結果與精確解對比

(9)

圖5給出了激波作用R22氣泡后典型過程的計算結果(密度云圖)與試驗結果對比，R22氣泡從開始壓扁，到繼續(xù)變形以及形成上下兩個典型漩渦的過程。計算結果與試驗結果吻合較好，證明了本文的計算模型能夠較為精確地捕捉激波的傳播以及多相流體界面的運動過程。

圖5 二維氣泡運動過程的密度云圖計算與試驗對比

5 水下近壁面爆炸氣泡射流現(xiàn)象

裝藥爆轟采用瞬時爆轟模型，流場初始條件為

(10)

可以看出，氣泡在21.2 ms時刻附近膨脹到最大體積，接下來氣泡開始收縮過程，周圍流體對氣泡進行壓縮。由37.1 ms時刻流場的壓力云圖可知，在氣泡收縮過程中，左側氣泡壁與固壁之間流場壓力較低，兩者之間的相對運動處于較為平衡狀態(tài)，而右側氣泡壁附近流場中存在局部高壓區(qū)，使得氣泡壁面加速向左運動。由于氣泡表面運動的非均衡性，在氣泡收縮過程中右側氣泡壁逐漸出現(xiàn)凹陷，開始產生射流形態(tài)，并沖擊左側氣泡壁。在39.3 ms附近射流開始沖擊左側氣泡壁，出現(xiàn)水錘現(xiàn)象，在壁面附近產生局部高壓，開始向外輻射壓力。之后爆炸氣泡也由開始的單連通域演變成多連通域，出現(xiàn)渦環(huán)現(xiàn)象，接下來開始氣泡的第二次脈動過程。

圖6 近壁面水下爆炸計算模型

圖7 近壁面水下爆炸的密度與壓力演化過程

圖8所示為在水射流出現(xiàn)階段爆炸氣泡的等效半徑與壁面中心點處流場的壓力時程曲線，圖中標出了與圖7對應的部分典型時刻在壓力曲線上的位置。從圖9可以看出，39.3 ms時刻射流開始沖擊左側氣泡壁時，水射流壓力仍處于上升階段，氣泡仍處于壓縮階段。39.8 ms水射流壓力達到峰值壓力46.8 MPa，此時錐型水柱主體部分作用到左側氣泡壁。之后隨著水射流能量的逐步消耗，壓力也逐漸下降。在壓力下降過程中出現(xiàn)的部分震蕩現(xiàn)象，結合圖7可以看出，此時氣泡界面并不完整，出現(xiàn)多連通域。在射流撞擊產生沖擊波及其傳播的過程中，容易疊加部分高頻震蕩。圖9為39.3 ms時刻流場中水平方向的速度分布云圖，最大速度超過了280 m/s。

圖8 爆炸氣泡等效半徑與壁面中心點壓力的曲線

圖9 流場中水平方向速度分布云圖

6 結 論

本文基于多相可壓縮流體的 five -equation 計算模型，引入界面函數(shù)壓縮和密度修正等技術來克服水射流演化過程中出現(xiàn)的界面彌散效應，獲得了近壁面水下爆炸爆炸氣泡運動的典型過程，獲得了水射流沖擊壁面產生的直接載荷特性。研究結果表明，水射流載荷的峰值壓力較大，且有效作用時間較長，其對于附近結構的沖擊損傷不能忽略。雖然目前針對水下爆炸氣泡潰滅及其壁面效應的實驗結果較為豐富，但由于水射流方向的隨機性以及測試技術的限制，目前對于水射流直接載荷的測量數(shù)據(jù)少見。后續(xù)擬在此方面開展相關研究，進一步考核驗證本文的計算模型，使之更加精確適用。