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

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

1 引 言

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

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

2 計(jì)算模型

2.1 控制方程

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

(1)

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

2.2 狀態(tài)方程

為使系統(tǒng)(1)封閉，需要補(bǔ)充流體的狀態(tài)方程。在流體計(jì)算模型中普遍采用的是剛性氣體方程(SG)，能夠描述多種液相或者氣相流體的運(yùn)動(dòng)特性[16,17]。近年來(lái)NASG方程逐漸得到采納，該方程修正了SG方程對(duì)液相介質(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é)算例涉及到水中空化問(wèn)題所采用的計(jì)算參數(shù)均采用表1的參數(shù)[14]。

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

2.3 空化模型

控制方程(1)并未考慮各相之間發(fā)生的相變過(guò)程，是從純物理幾何運(yùn)動(dòng)的守恒角度來(lái)描述流體的守恒運(yùn)動(dòng)過(guò)程。對(duì)于流體中發(fā)生液相與其對(duì)應(yīng)蒸汽相介質(zhì)之間物質(zhì)與能量轉(zhuǎn)換的情況，可以參考化學(xué)反應(yīng)過(guò)程中判斷過(guò)程進(jìn)行方向的判據(jù)，認(rèn)為在達(dá)到平衡態(tài)之后液相與蒸汽相之間的吉布斯(Gibbs)自由能相等。同時(shí)結(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ù)值方法簡(jiǎn)介

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

4.1 含氣率較高的激波管問(wèn)題測(cè)試

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

圖2 T=1 ms時(shí)刻流場(chǎng)密度、壓力、速度、溫度、液相與氣相質(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èn)題測(cè)試

計(jì)算域、網(wǎng)格數(shù)、初始間斷位置以及兩側(cè)初始?jí)毫ν?.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時(shí)刻各物理量的分布?？梢钥闯?，即使在水蒸氣含量極低的條件下，計(jì)算模型也能對(duì)右側(cè)沖擊波傳播過(guò)程中引起的蒸汽冷凝以及左側(cè)稀疏波傳播過(guò)程中引起液態(tài)蒸發(fā)等現(xiàn)象進(jìn)行較好的捕捉。

4.3 雙稀疏波問(wèn)題測(cè)試

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

圖3 T=1.5 ms時(shí)刻流場(chǎng)密度、壓力、速度、溫度、液相與氣相質(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時(shí)刻流場(chǎng)密度、壓力、速度、溫度、液相與氣相質(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)象的模擬

采用二維軸對(duì)稱模型研究近水面爆炸片空化現(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。計(jì)算域取8 m×14 m,其中水深11m，空氣厚度3 m，采用均勻網(wǎng)格，步長(zhǎng)0.02 m。圖5為t=1.5 ms， 3.2 ms 和4.6 ms時(shí)刻流場(chǎng)中密度、水蒸氣質(zhì)量分?jǐn)?shù)、溫度以及壓力的分布?？梢钥闯?，從初始的沖擊波傳播到水面后反射稀疏波，到爆炸氣體與水面之間出現(xiàn)低壓區(qū)域，以及空化閉合等三個(gè)典型時(shí)刻。圖5中每幅云圖的觀察窗口均為[-8,8]×[1.5,14] m2。

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

圖5 T=1.5 ms，3.2 ms和4.6 ms時(shí)刻流場(chǎng)密度、蒸汽相質(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 對(duì)稱軸上距水面不同間距處的壓力時(shí)程曲線Fig.6 Pressure time -history curves for different vertical distance from the initial water surface on the axisymmetric axis

5 結(jié) 論

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