基于相變松弛法則的水下爆炸空化模型研究

余 俊， 王海坤， 劉建湖， 汪 俊， 毛海斌， 周章濤

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

1 引 言

空化作為水下爆炸過程中的一種重要現(xiàn)象，其對水中的沖擊波傳播、爆炸氣泡運(yùn)動以及結(jié)構(gòu)物的沖擊損傷都會產(chǎn)生重要影響。對于水下爆炸過程中空化現(xiàn)象的研究，早期可以追溯到Taylor[1]針對水下爆炸流固耦合過程提出的Taylor平板模型，對空化的產(chǎn)生機(jī)理進(jìn)行了理論探索。此后，文獻(xiàn)[2-5]基于各種假設(shè)對水下爆炸空化現(xiàn)象進(jìn)行了理論分析與解釋，并建立了部分經(jīng)驗(yàn)公式。由于水下爆炸空化過程的復(fù)雜性，理論模型研究對象基本都是針對簡單理想工況，其應(yīng)用范圍受到了較大的限制。在試驗(yàn)研究方面，文獻(xiàn)[6-8]對近水面爆炸空化現(xiàn)象進(jìn)行了記錄與觀察，為水下爆炸空化的機(jī)理研究提供了部分寶貴的試驗(yàn)數(shù)據(jù)。理論研究的局限性以及實(shí)驗(yàn)研究的不可預(yù)見性，使得水下爆炸空化的數(shù)值模擬研究顯示出較大的優(yōu)勢。目前空化模型主要分為 one -fluid模型與 two -fluid模型兩大類。one -fluid模型主要以cut-off模型[9]、isentropic模型[10]和Schmidt模型[11]等為典型代表，其主要思想是當(dāng)液體中的壓力小于飽和壓力，認(rèn)為就會出現(xiàn)空化，其中各種模型的具體差別體現(xiàn)在空化流場中壓力或者密度計算的復(fù)雜度不同。two -fluid模型則主要以文獻(xiàn)[12-15]提出的4-euqation和6-equation等模型為主，將可壓縮流體控制方程與空化相變模型相結(jié)合，通過追蹤空化質(zhì)量分?jǐn)?shù)的演化來捕捉空化的產(chǎn)生、發(fā)展與潰滅過程。Chiapolino等[12]提出了汽-液兩相平衡態(tài)計算的簡化方法，避免了繁雜的迭代計算過程，提出了空化相變計算效率。Zhang[19]將five-equation與相變模型結(jié)合來模擬空化。

受到認(rèn)識水平以及計算方法的限制，目前水下爆炸空化計算領(lǐng)域使用的基本均是 one -fluid模型。與 one -fluid模型相比，two -fluid模型要復(fù)雜得多，需要充分考慮空化熱力學(xué)運(yùn)動過程中的質(zhì)量與熱量轉(zhuǎn)換，對于計算的守恒性要求更嚴(yán)格。但是 two -fluid模型能夠持續(xù)地追蹤蒸汽相的演化過程，從邏輯上來說該方法更加科學(xué)合理。因此，本文擬采用4-equation模型來描述多相流運(yùn)動，其中汽-液兩相之間的相變過程引入松弛法則進(jìn)行處理。在此基礎(chǔ)上對上述模型進(jìn)行了數(shù)值測試，最后將其應(yīng)用到近水面爆炸空化研究領(lǐng)域。

2 計算模型

2.1 控制方程

對于水下爆炸運(yùn)動等瞬態(tài)響應(yīng)過程，可以忽略傳熱、粘性以及化學(xué)反應(yīng)等影響，采用基于等壓和等速度假設(shè)條件下的多相可壓縮流體計算模型4-equation 來描述[12]：

(1)

控制方程中參數(shù)k表示混合流體中各相成分，為了計算模型表達(dá)簡便，本文約定，k=1表示液相成分，k=2表示蒸汽相成分，k=3,…,N表示其他氣相或液相成分，但其中的氣相不會轉(zhuǎn)換為其對應(yīng)的液相，液相也不會轉(zhuǎn)換為其對應(yīng)的氣相。統(tǒng)一令N=3，且k=3只表示氣相，對于其他情況可以簡單類推。

2.2 狀態(tài)方程

為使系統(tǒng)(1)封閉，需要補(bǔ)充流體的狀態(tài)方程。在流體計算模型中普遍采用的是剛性氣體方程(SG)，能夠描述多種液相或者氣相流體的運(yùn)動特性[16,17]。近年來NASG方程逐漸得到采納，該方程修正了SG方程對液相介質(zhì)中分子排斥力效應(yīng)描述的不足[18]，可表示為

(2)

(3)

式中 下標(biāo)sat表示飽和態(tài)，參數(shù)A,B,C,D和E分別為

(4)

式中Cp,k為k相流體定壓比熱容。液態(tài)水、水蒸汽以及空氣的NASG狀態(tài)方程參數(shù)列入表1，第4節(jié)算例涉及到水中空化問題所采用的計算參數(shù)均采用表1的參數(shù)[14]。

表1 液態(tài)水、水蒸氣和空氣的NASG方程參數(shù) (國際單位)Tab.1 NASG coefficients for liquid water,vapor water and air

2.3 空化模型

控制方程(1)并未考慮各相之間發(fā)生的相變過程，是從純物理幾何運(yùn)動的守恒角度來描述流體的守恒運(yùn)動過程。對于流體中發(fā)生液相與其對應(yīng)蒸汽相介質(zhì)之間物質(zhì)與能量轉(zhuǎn)換的情況，可以參考化學(xué)反應(yīng)過程中判斷過程進(jìn)行方向的判據(jù)，認(rèn)為在達(dá)到平衡態(tài)之后液相與蒸汽相之間的吉布斯(Gibbs)自由能相等。同時結(jié)合控制方程(1)的假定，可以認(rèn)為系統(tǒng)達(dá)到平衡態(tài)后滿足關(guān)系式[19]為

(5)

將式(2,3)代入式(5)，可轉(zhuǎn)換為

(6)

式中ppartial為氣體混合物中蒸汽相的分壓，由式(6)可知其與蒸汽相的摩爾質(zhì)量分?jǐn)?shù)成正比。非線性方程(6)可以采用迭代方法進(jìn)行求解，如Newton-Raphson方法。

3 數(shù)值方法簡介

控制方程(1)與狀態(tài)方程(2,3)以及空化相變方程(6)共同構(gòu)成了帶有相變轉(zhuǎn)換的多相可壓縮流體系統(tǒng)方程組。該類計算模型的數(shù)值離散可以利用相變松弛法則進(jìn)行分步求解[19]。第一步求解齊次雙曲型方程組(1)，需要結(jié)合狀態(tài)方程(2)，利用二階MUSCL-Hancock方法以及HLLC近似黎曼求解器。通過該步獲得中間步變量(υ,e,p,T,Yk)。第二步采用Newton-Raphson迭代方法求得相變轉(zhuǎn)換平衡態(tài)時的狀態(tài)量(p*,T*,Yk =1,2)。

(7)

式中 參數(shù)ξi和ξ定義如下。

(8)

圖1 常溫常壓下氣液兩相混合流體聲速隨氣體體積分?jǐn)?shù)的變化Fig.1 Sound speed of gas-liquid two -phase mixture with gas volume fraction at room temperature and pressure

4 數(shù)值測試及算例

4.1 含氣率較高的激波管問題測試

將整個流場中的蒸汽相初始質(zhì)量分?jǐn)?shù)設(shè)置為0.2，液相與氣相的質(zhì)量分?jǐn)?shù)分別為0.1和0.7。計算域?yàn)閇0,1 m]，均勻劃分為400個網(wǎng)格。流場中密度和壓力的初始間斷面位置均為x=0.5 m，左右側(cè)壓力分別為2 bar和1 bar，兩側(cè)密度分別為1.94 kg/m3和1.02 kg/m3。圖2為t=1 ms時刻密度、壓力、速度、溫度、液相與蒸汽相質(zhì)量分?jǐn)?shù)分布曲線，其中P -T和no P -T曲線分別表示有/無相變轉(zhuǎn)換的計算工況，initial表示初始條件。從圖2密度與壓力分布可以看出，考慮相變轉(zhuǎn)換后，兩側(cè)的沖擊波與稀疏波波速均有所增加，溫度與質(zhì)量分?jǐn)?shù)分布差異比較明顯?？梢园l(fā)現(xiàn)向右側(cè)傳播的沖擊波壓縮低壓流體后，流體中的蒸汽相質(zhì)量分?jǐn)?shù)有所增加，說明混合流體中的液相發(fā)生了相變，部分發(fā)生蒸發(fā)。同時向左側(cè)傳播的稀疏波使得左側(cè)高壓流體中水蒸氣發(fā)生冷凝，液相質(zhì)量分?jǐn)?shù)有所增加。

圖2 T=1 ms時刻流場密度、壓力、速度、溫度、液相與氣相質(zhì)量分?jǐn)?shù)分布Fig.2 Distribution diagram of density,pressure,velocity,temperature,liquid phase and vapor phase mass fraction of flow field at T=1 ms

4.2 含氣率極小的激波管問題測試

計算域、網(wǎng)格數(shù)、初始間斷位置以及兩側(cè)初始壓力同4.1節(jié)。左側(cè)蒸汽相與空氣的質(zhì)量分?jǐn)?shù)分別為9.45×10-8與10-5，右側(cè)分別為1.92×10-7與10-5，兩側(cè)密度分別為1051.4 kg/m3與1046.5 kg/m3。圖3為1.5 ms時刻各物理量的分布?？梢钥闯觯词乖谒魵夂繕O低的條件下，計算模型也能對右側(cè)沖擊波傳播過程中引起的蒸汽冷凝以及左側(cè)稀疏波傳播過程中引起液態(tài)蒸發(fā)等現(xiàn)象進(jìn)行較好的捕捉。

4.3 雙稀疏波問題測試

計算域、網(wǎng)格數(shù)、初始間斷位置以及兩側(cè)初始壓力同4.1節(jié)。初始流場中液相、蒸汽相和氣相所占質(zhì)量分?jǐn)?shù)均為 0.0304， 0.5 和 0.4696，壓力為 105Pa。初始間斷位于x=0.5 m處，間斷面左側(cè)速度和右側(cè)速度分別為-20 m/s和20 m/s。計算網(wǎng)格N=400，圖4為t=0.5 ms時刻流場中密度、壓力、速度、溫度、液相和蒸汽相質(zhì)量分?jǐn)?shù)的分布?？梢钥闯觯诳紤]了相變過程之后，流場的密度和壓力等分布變化較大，主要還是由于其中的氣相與液相之間發(fā)生了相變轉(zhuǎn)換，兩相質(zhì)量分?jǐn)?shù)都會產(chǎn)生變化，從而引起各物理量的重新分布。

圖3 T=1.5 ms時刻流場密度、壓力、速度、溫度、液相與氣相質(zhì)量分?jǐn)?shù)分布Fig.3 Distribution diagram of density,pressure,velocity,temperature,liquid phase and vapor phase mass fraction of flow field at T=1.5 ms

圖4 T=0.5 ms時刻流場密度、壓力、速度、溫度、液相與氣相質(zhì)量分?jǐn)?shù)分布Fig.4 Distribution diagram of density,pressure,velocity,temperature,liquid phase and vapor phase mass fraction of flow field at T=0.5 ms

4.4 近水面爆炸空化現(xiàn)象的模擬

采用二維軸對稱模型研究近水面爆炸片空化現(xiàn)象。其中爆炸氣體半徑0.5 m，位于水下 3 m 位置，初始密度和壓力分別為1343 kg/m3和4.29×108Pa，內(nèi)部的水蒸氣和氣體質(zhì)量分?jǐn)?shù)分別為0.05和0.95。水面上方空氣的初始密度和壓力分別為1.3 kg/m3和105Pa，內(nèi)部的水蒸氣和氣體的質(zhì)量分?jǐn)?shù)分別為0.001和0.999。水的初始密度和壓力分別為1050 kg/m3和105Pa，內(nèi)部的水蒸氣和氣體的質(zhì)量分?jǐn)?shù)分別為10-7和4.73×10-6。計算域取8 m×14 m,其中水深11m，空氣厚度3 m，采用均勻網(wǎng)格，步長0.02 m。圖5為t=1.5 ms， 3.2 ms 和4.6 ms時刻流場中密度、水蒸氣質(zhì)量分?jǐn)?shù)、溫度以及壓力的分布?？梢钥闯觯瑥某跏嫉臎_擊波傳播到水面后反射稀疏波，到爆炸氣體與水面之間出現(xiàn)低壓區(qū)域，以及空化閉合等三個典型時刻。圖5中每幅云圖的觀察窗口均為[-8,8]×[1.5,14] m2。

圖6為縱向?qū)ΨQ軸上距水面分別為0.2 m,0.4 m 和0.6 m處的水中壓力時程曲線。可以看出，在2 ms～4 ms期間，空化處于產(chǎn)生與發(fā)展階段，4 ms之后空化開始受到壓縮，逐漸潰滅，最大潰滅壓力達(dá)到7.6 MPa，約為沖擊波峰值壓力的15%。另外空化潰滅載荷的有效作用時間在1 ms～2 ms，是沖擊波有效作用時間的2倍以上，可見其對結(jié)構(gòu)的沖擊效應(yīng)不能忽略。

圖5 T=1.5 ms，3.2 ms和4.6 ms時刻流場密度、蒸汽相質(zhì)量分?jǐn)?shù)、溫度和壓力云圖Fig.5 Pseudo -color plots of the mixture density,vapor mass fraction,temperature and pressure at different times T =1.5 ms,3.2 ms and 4.6 ms

圖6 對稱軸上距水面不同間距處的壓力時程曲線Fig.6 Pressure time -history curves for different vertical distance from the initial water surface on the axisymmetric axis

5 結(jié) 論

為了捕捉水下近水面爆炸過程中的空化現(xiàn)象，本文在多相可壓縮流體模型基礎(chǔ)上，考慮汽-液兩相流體亞平衡狀態(tài)下質(zhì)量與熱量交換過程中的熱力學(xué)-化學(xué)平衡機(jī)制，實(shí)現(xiàn)了對空化相變的自動捕捉。控制方程采用分步法求解，首先利用 MUSCL-Hancock與HLLC黎曼求解器求解齊次雙曲型方程，然后再利用迭代方法求解相變松弛方程。對該計算模型進(jìn)行了多種工況下的空化相變測試，對于汽-液兩相之間的蒸發(fā)與冷凝等相變過程具有較好的捕捉能力。最后將該模型應(yīng)用于水下近水面爆炸空化現(xiàn)象的模擬，展示了水下爆炸空化的整個過程，捕捉了完整的空化閉合壓力，本文計算方法與研究結(jié)果對提高水下爆炸空化機(jī)理的認(rèn)識及其深入研究具有借鑒和指導(dǎo)價值。