余 俊，王海坤，劉國振，郝 軼，張得揚

（1.中國船舶科學研究中心，江蘇無錫 214082；2.深海技術科學太湖實驗室，江蘇無錫 214082）

0 引 言

隨著現(xiàn)代海戰(zhàn)武器的高速發(fā)展，實戰(zhàn)條件下艦船遭受的威脅越來越大，對艦船水下沖擊與防護性能的研究受到了高度的重視。水面艦艇除了會遭受船底爆炸沖擊損傷之外，近水面爆炸也會產生重要的威脅[1]。近水面爆炸現(xiàn)象的本質特征屬于多相流體（氣-液-汽）相互作用及其與結構的耦合運動過程，其非線性耦合效應主要體現(xiàn)在兩個重要方面：首先是爆轟氣體產物、水、空氣等多組分流體之間的對流運動,其可歸納為多相可壓縮流體界面捕捉問題；其次是自由面的存在使得流場中會出現(xiàn)稀疏波的傳播，水中會出現(xiàn)典型空化現(xiàn)象,可以歸結為相變問題。要想準確理解和掌握近水面爆炸過程中的多流體非線性耦合效應，需要同時重點處理好多相可壓縮流體界面捕捉和相變問題。

目前國內外已經對近水面爆炸現(xiàn)象進行了大量的研究[2-4]。在實驗研究方面，Marcus等[5]開展了水下近自由面爆炸試驗研究，發(fā)現(xiàn)水下爆炸產生的片空化區(qū)域在閉合潰滅時會產生較大的壓力，獲得了近水面爆炸過程中包括空化效應的沖擊波和氣泡階段的壓力曲線;Brett等[6]開展了一系列水下圓柱殼近距離爆炸試驗，對沖擊波和空化聯(lián)合加載下圓柱殼的變形進行了測量分析，發(fā)現(xiàn)空化載荷的幅值與沖擊波幅值相當，而且空化潰滅載荷造成的殼體變形要超過整個沖擊波階段產生的總變形；Kleine等[7]利用紋影法拍攝到了近水面爆炸試驗過程中液體空化的密度變化過程圖象，較為清晰地展示了爆炸氣泡膨脹以及空化域內部演化圖像；Cui 等[8]開展了一系列不同邊界條件下近水面小藥量藥包爆炸試驗，獲得了爆炸氣泡與自由面運動的完整過程。由于實驗研究獲得的數(shù)據(jù)有限，在穩(wěn)定性和重復性上難以滿足要求，難以獲得流場的精細特征。數(shù)值模擬能夠較為全面地展示流場的精細特征，使得其在近水面爆炸研究中具有一定的優(yōu)勢。這里著重介紹基于多相可壓縮流體計算模擬的研究情況，對于基于聲學單元和邊界單元的數(shù)值方法不作討論。自從Aanhold 等[9]提出cut-off 空化模型之后，由于其使用簡單方便，被廣泛應用在水下近自由面或近結構爆炸分析中。Wardlaw等[10]利用cut-off模型分析了充水圓筒中心藥包爆炸過程中沖擊波在圓筒內部產生的局部空化現(xiàn)象，獲得了空化潰滅載荷的典型特征；Shukla[11]、Yu[12]等采用cut-off 空化模型與五方程相結合的方法對帶有空化效應的近水面爆炸多流體運動問題進行了研究；Wu 等[13]采用cut-off模型與流體局部間斷迦流金（LDG）法相結合分析了水下遠場沖擊波與自由面作用產生空化的過程。為了克服cut-off模型無法描述空化域內液相和汽相組分的問題，Liu[14]、Schmidt[15]和Xie[16]等相繼提出了isentropic 模型、Schmidt 模型和Modified Schmidt模型。這些模型本質上都屬于單流體空化模型，在流體壓力降低到飽和態(tài)壓力之前均視為純流體，當降到飽和態(tài)壓力之后，各種模型的具體差別體現(xiàn)在空化域內汽液混合流體狀態(tài)方程的構造方法不同。除了上述單流體空化模型之外，雙流體模型也在不斷發(fā)展，主要以Chiapolino[17-18]、Saurel[19]、Pelanti[20-21]等提出的四方程、六方程等模型為主，該模型通過持續(xù)追蹤空化質量分數(shù)的演化來捕捉空化的產生、發(fā)展與潰滅過程。雙流體模型則將初始流體視為液相和汽相的某種混合，在流場發(fā)生改變時（如沖擊波傳播和氣泡運動等包含間斷時）液-汽兩相之間會發(fā)生對流運動，同時還伴隨有質量和熱量交換的相變過程，空化被視為流體區(qū)域內部蒸汽含量不斷累積的過程。因此雙流體模型描述的空化演化過程更符合物理本質和思維邏輯，但雙流體模型相比單流體模型，雙流體計算過程非常復雜繁瑣，導致目前水下爆炸領域內絕大部分的空化研究都采用單流體模型。

本文擬采用Chiapolino 等[18]提出的基于多相可壓縮流體四方程的雙流體空化模型來重點研究近水面爆炸過程中的多相流運動。首先對計算模型的控制方程及其數(shù)值方法進行了簡要介紹，然后采用一維激波管問題對計算模型中的相變轉換以及激波捕捉功能進行了考核，最后針對近水面爆炸現(xiàn)象，將計算結果與試驗結果進行對比，并重點分析空化域內部的典型物理特性。

1 計算模型簡介

1.1 控制方程

不考慮粘性和熱傳導效應的多相可壓縮流體運動方程可表示為[18]

式中，ρ、u、p和E分別代表混合流體的密度、速度矢量、壓力和單位質量總能，E=e+0.5u·u，這里e為單位質量內能；Yk代表混合流體中各相介質的質量分數(shù)?？刂品匠蹋?）中參數(shù)k代表了混合流體中各相成分，為了計算模型中表達簡便起見，這里約定如下：k=1代表液相成分；k=2代表蒸汽相成分；k=3,…，N代表其它不發(fā)生相變的氣相或液相成分。本文中令k=3 代表空氣（不可冷凝），k=4 代表爆轟氣體產物。

在流體計算模型中普遍采用的是SG（Stiffed-Gas）狀態(tài)方程，該方程能夠描述多種液相或者氣相流體的運動特性。近年來NASG（Nobel-Able Stiffed-Gas）狀態(tài)方程逐漸被采納，該方程修正了SG方程對液相介質中分子排斥力效應描述的不足[22]。該方程形式為

式中，?k=1/ρk、Tk、gk、ck分別代表k相流體介質的比體積、溫度、吉布斯（Gibbs）自由能與聲速，參數(shù)γk、、Ck、qk，qk'、bk等可以通過流體熱動力學試驗中與溫度相關的飽和態(tài)壓力、各相比體積、各相的焓等數(shù)據(jù)來擬合得到[22]。在飽和平衡態(tài)下液相與汽相的Gibbs 自由能相等，即g1=g2，化簡可得飽和態(tài)下溫度與壓力的關系式：

式中下標“sat”代表飽和態(tài)，參數(shù)As、Bs、Cs、Ds、Es分別為

式中Cp,k表示第k相流體定壓比熱容。液態(tài)水、水蒸汽以及空氣的NASG 狀態(tài)方程參數(shù)如表1所示，其中W為物質的摩爾質量。本文算例涉及到水中空化問題所采用的計算參數(shù)均采用表1中的參數(shù)[22]。

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

1.2 相變轉換模型

控制方程（1）并未考慮到各相之間發(fā)生的相變過程，對于流體中發(fā)生液相與其對應汽相之間的物質與熱量轉換，可以參考化學反應過程中判斷反應方向的判據(jù)，認為在達到平衡態(tài)之后液相及其蒸汽相之間的吉布斯自由能相等。同時結合控制方程（1）的假定，認為系統(tǒng)達到平衡態(tài)后滿足如下關系式[18]：

將公式（2）和（3）代入，得

這里ppartial代表氣體混合物中蒸汽相的分壓，由公式可知其與蒸汽相的摩爾質量分數(shù)成正比[18]。非線性方程（6）～（8）可以采用迭代方法進行求解，如Newton-Raphson方法。

1.3 數(shù)值離散

控制方程（1）與狀態(tài)方程（2）以及空化相變方程（6）～（8）共同構成了考慮相變效應的多相可壓縮流體控制方程組，可以利用相變松弛法則進行分步求解[23]。第一步求解齊次雙曲型方程組（1），需要結合狀態(tài)方程（2），利用二階MUSCL-Hancock方法以及HLLC近似黎曼求解器。通過該步獲得中間變量（v,e,p,T,Yk）；第二步采用Newton-Raphson 迭代方法求得相變轉換平衡態(tài)時的狀態(tài)量（p*,T*,Yk=1,2）。

2 計算與討論

2.1 一維激波管測試與驗證

本算例主要用來與Chiapolino 等相同計算工況[18]進行對比驗證。其中整個流場的蒸汽相初始質量分數(shù)設置為0.2，液相與氣相的質量分數(shù)分別為0.1和0.7，計算域為[0,1]m。流場中密度和壓力的初始間斷面位于x=0.5 m，左、右側壓力分別為2×105Pa和105Pa，兩側密度分別為1.94 kg/m3、1.02 kg/m3。

圖1 為t=1 ms 時刻流場中密度、壓力、速度、溫度、液相與蒸汽相質量分數(shù)分布曲線。其中P-T 和no P-T 曲線分別代表是否考慮相變轉換時的計算工況，initial代表初始條件，Chiapolino 代表對比的數(shù)據(jù)（網格數(shù)為100）。同時為了考察計算模型的網格收斂性，分別采用100、200 和400 個網格單元均勻劃分計算域。由密度和壓力分布曲線可知，考慮液相和汽相之間的相變轉換之后，向右傳播的壓縮波和向左傳播的稀疏波的波速都稍有下降。由溫度曲線可知，向右傳播的壓縮波波后溫度在兩種工況下都有所升高，但波后穩(wěn)定區(qū)的溫度在考慮相變轉換時為346.3K,而不考慮相變時為362.6 K，差異較大。同理，向左傳播的稀疏波波后穩(wěn)定區(qū)的溫度都有所下降，但是考慮相變轉換后為344.7 K，而不考慮時則為329.4 K。這是由于在壓縮波傳播的波陣面以及稀疏波傳播過程中液相與汽相發(fā)生了質量和熱量交換，在達到平衡態(tài)的過程中溫度存在變換。而汽-液兩相的質量轉換可以參考質量分數(shù)曲線得知，壓縮波過后流體中液相發(fā)生了蒸發(fā)，使得汽相質量分數(shù)增加。該計算結果與通常理解的沖擊波作用后由于壓力升高容易引起液相發(fā)生冷凝的認識相矛盾，這里需要結合溫度變化進行聯(lián)合分析。由于考慮了相變轉換，沖擊波作用后波后溫度由初始的337.5 K 上升到346.3 K,雖然波后壓力也由初始的105Pa上升到1.4×105Pa，但是從飽和態(tài)曲線P-T（公式（3））來分析會發(fā)現(xiàn)溫度升高引起的蒸發(fā)效應要大于壓力升高帶來的冷凝效應，最終效果就是壓縮波后處于蒸發(fā)狀態(tài)。同理，向左傳播稀疏波的最終效果是波后處于冷凝狀態(tài)。由圖1 可知，隨著網格數(shù)的增加，計算結果的收斂性較好，本文的計算結果與Chiapolino的結果吻合較好。

圖1 1 ms時刻流場各物理量分布曲線圖Fig.1 Diagram of different variables in fluid at t=1 ms

2.2 水下近自由面爆炸現(xiàn)象

本節(jié)利用上面介紹的計算模型與Cui 等[8]的試驗觀察結果進行對比。該試驗在邊長2 m 的立方體水槽內開展，采用等效于5.2 g TNT 的藥包起爆，藥包中心位于水下0.13 m 處。藥包采用等效爆轟模0.0型09近m似，爆處轟理氣，根體據(jù)產G物e內e r-部H密u n度ter和理壓論力模分型別[24]以為及16數(shù)06值kg經/m驗3和，初10始9P球a，形采爆用轟理氣想體氣產體物狀半態(tài)徑方設程置，γ為=1R.80=，CV=695，W=290?？諝?、水和爆炸氣泡內的初始蒸汽質量分數(shù)分別為10-9、10-8和10-6?？諝狻⑺捅馀莸某跏紲囟确謩e設置為295 K、295 K 和1120 K,其中空氣和爆轟氣體產物處于過熱狀態(tài)，水處于飽和狀態(tài)。空氣、水和爆炸氣泡內的初始蒸汽體積分數(shù)分別為1.56×10-9、1.39×10-5和8.07×10-7。為了節(jié)省計算資源，采用二維軸對稱模型進行模擬，如圖2 所示。其中y軸為對稱旋轉軸，y=0 為初始水平面位置，藥包位于(0,-0.13)m處。計算域采用非均勻結構性網格劃分，如圖2 所示（圖中顯示的網格密度為實際的1/2）。其中區(qū)域[0, 0.4]×[-0.6, 0.2] m2內采用均勻網格，dx=dy=0.001 m，區(qū)域外采用等比數(shù)列格式的遞增步長網格，等比為1.02。整個計算域網格為650×1200，計算域大小為[0,7.34]×[-7.5,1.1]m2。y軸邊界設置為對稱條件，其它如空氣層上方、水域下方以及右側邊界均設置為無反射條件，CFL設置為0.5。

圖2 近水面爆炸二維軸對稱計算模型與網格劃分示意圖Fig.2 Schematic of 2D axisymmetric model and mesh size of underwater explosion near free surface

圖3 顯示了在考慮相變條件下t=0.084 ms、0.168 ms、0.334 ms、0.5 ms、0.666 ms、1.0 ms、4.176 ms 和6.176 ms 時刻流場中的密度、壓力和蒸汽體積分數(shù)分布云圖，其中所有圖片采用統(tǒng)一大小的觀察窗，尺寸為[-0.6,0.6]×[-0.7,0.35]m2。由于圖3 中的各物理量的內部分布跨度非常大，未單獨列出各云圖的刻度線，壓力云圖中的黑色虛線表示對應時刻爆炸氣泡界面的位置。0.084 ms 時刻，初始沖擊波剛到達自由面不久，其產生的反射稀疏波使得自由面下方附近出現(xiàn)蒸汽聚集，圖中對應時刻蒸汽體積分數(shù)云圖顯示的局部高亮區(qū)域內數(shù)值范圍為0.4‰～2.5‰。隨著稀疏波不斷向下傳播，0.168 ms 時刻稀疏波已經達到爆炸氣泡上壁面，并反射壓縮波，可以發(fā)現(xiàn)此時爆炸氣泡表面附近的局部范圍處于相對高壓區(qū)。隨著稀疏波在爆炸氣泡表面反射成壓縮波以及自由面不斷反射稀疏波，旋轉軸上的蒸汽含量不斷降低（對應空泡潰滅過程），而旋轉軸外側的蒸汽含量不斷增加（對應空化產生過程），如圖3中的0.334 ms和0.500 ms所示。隨著空泡的不斷產生和潰滅，空泡中心區(qū)域逐漸變薄，并向外擴張，如圖3中的0.666 ms和1.0 ms所示。在4.167 ms和6.167 ms時刻空泡基本上完全消失，但爆炸氣泡一直處于向外膨脹過程，并不斷抬高水面位置。在整個過程中空化由開始的單連通域不斷演化為多連通域，顯示出渦環(huán)形狀，渦環(huán)逐漸向外擴展，并逐漸變薄直至最終完全消失。

圖3 典型時刻流場的密度、壓力和蒸汽體積分數(shù)分布云圖Fig.3 Contour plots of the density,pressure and vapor volume fractions at typical times

在上述整個時間段內，自由面上方附近的空氣中蒸汽含量一直比較高，并且包含蒸汽的區(qū)域不斷擴展。這是由于水中初始沖擊波作用自由面后向空氣透射壓縮波、水中的液相和汽相同時向空氣中擴散的結果，在擴散過程中仍然存在液-汽兩相之間的相變轉換。這是因為水面上方初始空氣中是不包含液相水，而且水蒸汽含量極低（10-9），而圖3中蒸汽體積分數(shù)云圖中紅色區(qū)域代表自由面上方，其內部的液相水的質量分數(shù)范圍在0.6～0.8，水蒸汽質量分數(shù)范圍在0.1～0.3，剩下的為空氣。

文獻[8]中只記錄了通過觀察窗獲得的部分流場圖像，如圖4 中第一行所示。試驗獲得的圖片的中間縱剖面尺寸與數(shù)值模擬中的[-0.4,0.4]×[-0.68,0.2]m2基本相當。在第一行的試驗圖片中左下角兩位數(shù)值代表圖片序號，右下角為記錄的時刻（單位為ms），其中0.167 ms代表裝藥起爆時刻。由圖4可知，小藥量裝藥在近水面爆炸時在水面下方產生的空化域形態(tài)與“云空化”非常相似，水中的空泡體積非常小，形成云帶狀，其發(fā)展過程中能夠觀察到明顯的渦環(huán)形狀。為了與試驗進行對比，需要對“云空化”區(qū)域內的蒸汽含量進行判斷，選取蒸汽體積分數(shù)是一個比較合適的對象。這里假定當流體單元中的蒸汽體積分數(shù)大于0.5‰時就標記為空化單元，單元所在的區(qū)域即為空化域。有關空化域內蒸汽含量的判據(jù)目前仍然沒有統(tǒng)一的標準，這里選取的0.5‰只是數(shù)值計算上的一個參考。將圖3 中蒸汽體積分數(shù)云圖轉換為空化域紋影圖，如圖4 第二行所示。其中紅色區(qū)域為水下空化域，黃色區(qū)域為自由面上方的高含量蒸汽集中區(qū)域，綠色為空氣，藍色為水域，淺藍色為爆炸氣泡。經過對比可知，計算獲得的水中空化域與試驗觀察到的“云空化”范圍基本相當，兩者都清晰地展示了空化域形成渦環(huán)，并逐漸向外擴張和變薄的過程，最終完全潰滅。

圖4 典型時刻的試驗與計算結果對比Fig.4 Comparison of experimental and numerical results at typical times

圖5 顯示了空化域內的最大、最小和平均溫度變化以及空化域總體積時程曲線。由圖5 可知空泡總體積的膨脹過程與收縮過程的時間基本相當，整個空化運動周期約為1.14 ms?？栈蚺蛎浐褪湛s過程中內部最高溫度變化范圍為294.6～309.7 K，最低溫度變化范圍為283.7～294.6 K。膨脹到最大體積時空化域內部的最高和最低溫度分別為309.7 K 和283.7 K，但在空化域變化過程中其內部體積平均溫度基本恒定為294.6 K，與室溫相當。

圖6 為流場中三個坐標點處壓力時程曲線圖，分別為水下0.05 m 處距離爆心不同水平距離處的測點。由圖可知，三個測點處的空化潰滅過程中最大壓力分別為2.15 MPa、0.98 MPa 和0.92 MPa。圖7 為空化域內最大、最小和平均壓力變化時程曲線。由圖可知，空化域內最大壓力在體積膨脹階段的大部分時刻保持在0.018 MPa 左右，而在收縮階段則主要在0.018 MPa 上方振蕩，最大值達到0.029 MPa。空化域內最小壓力在中間的絕大部分時間內變化較為平緩，在3200～4800 Pa 范圍內變動，只是在空化剛開始產生和最后潰滅階段變化梯度較大?？栈騼润w積平均壓力變化范圍為9000～18 000 Pa，其變化趨勢與最小壓力變化趨勢基本保持一致。由此可見，空化域內壓力并不是維持在飽和態(tài)壓力不變，其具有分布不均勻性、變化梯度大等特點。同時結合圖5可知其內部溫度的變化范圍在283.7～309.7 K，空化域內部發(fā)生較為明顯的相變轉換現(xiàn)象。

圖5 空化域內最大、最小和平均溫度以及空化域體積的時程曲線Fig.5 Time history of the maximun,minimun,average temperatures and cavitation domain volumes in cavitation domain

圖6 三個測點處壓力時程曲線圖Fig.6 Time history of pressure at three measuring points

圖7 空化域內最大、最小和平均壓力時程曲線Fig.7 Time history of the maximum,minimum and average pressures in cavitation domain

圖8 為空化域內的蒸汽相體積分數(shù)的最大和平均值變化時程曲線，空化域內蒸汽相體積分數(shù)的最小值為前面設定參考值的0.5‰。由圖8 可知，空化域內的蒸汽相體積分數(shù)在0.5‰～173‰，其峰值并未出現(xiàn)在空化體積膨脹到最大的時刻。盡管空化整個運動過程中蒸汽相最大體積分數(shù)并不低，但是空化域內部整體體積平均值變化區(qū)間在0.5‰～10‰，說明蒸汽的平均含量比較低，這與前面介紹的實驗中觀察到的“云空化”較為一致，同時也說明了本文計算模型具有一定的合理性和精確性。圖9為空化域內密度的最大、最小和平均值變化時程曲線。由圖9可知，空化域內的最大密度變化范圍為1051～1061 kg/m3，體積平均密度變化范圍為1039～1051 kg/m3，最小密度變化范圍860～1061 kg/m3，其峰值也未出現(xiàn)在空化體積膨脹到最大的時刻。結合圖8可知，空化域內蒸汽相的含量總體上非常小，對于流體密度降低的程度有限，空化域內總體體積平均密度變化較小。

圖8 空化域內蒸汽體積分數(shù)?2最大值、平均值變化曲線Fig.8 Time history of the maximum and average vapor volume fractions ?2 in cavitation domain

圖9 空化域內密度的最大值、最小值和平均值的時程曲線Fig.9 Time history of the maximum,minimum and average densities in cavitation domain

3 結 論

針對近水面爆炸過程中出現(xiàn)的爆轟氣體產物、水以及空氣等多流體復雜運動，以及各種流體內部可能出現(xiàn)的液相水和汽相水之間的相變轉換等現(xiàn)象，本文采用了基于相變效應的多流體運動計算模型進行處理，獲得了如下結論：

（1）利用一維激波管問題對計算模型精確性和收斂性測試表明，本文的模型計算結果與Chiapolino的結果吻合較好，能夠精確地捕捉沖擊波與稀疏波的傳播過程，展示了沖擊波作用之后流體蒸發(fā)以及稀疏波作用之后流體冷凝的過程，最終的狀態(tài)往往是兩種過程的綜合作用的結果；

（2）近水面爆炸過程中爆炸氣泡運動過程的計算結果與試驗結果符合較好，說明本文建立的多流體耦合計算模型能夠捕捉多種界面之間的復雜運動過程。

（3）多流體計算模型中引入了相變轉換，可獲得近水面爆炸過程中水面下方的液體內部發(fā)生的空化現(xiàn)象；采用蒸汽體積分數(shù)0.5‰作為空化域識別的判據(jù)，數(shù)值模擬獲得的空化域范圍與試驗結果較為一致，空化域都是由開始的單連通域演化為渦環(huán)形態(tài)，并向外擴展直至完全潰滅。

（4）數(shù)值模擬結果發(fā)現(xiàn)，空化域在運動過程中其內部蒸汽的體積分數(shù)范圍為0.5‰～173‰，但內部平均值的范圍僅為0.5‰～10‰，說明蒸汽的平均體積含量比較低，這與實驗中觀察到的“云空化”現(xiàn)象具有一致性。同時該現(xiàn)象也可以從空化域內部密度的平均體積含量來得到印證，其變化范圍為1039～1051 kg/m3，說明空化域內部平均密度與常態(tài)下液相水的密度非常接近。另外，空化域內的壓力變化范圍為3200～29 000 Pa，溫度變化范圍為283.7～309.7 K，說明空化域內部的壓力和溫度一直處于動態(tài)變化之中，并不是維持在某個恒定的飽和態(tài)不變。

近水面爆炸過程中除了存在多流體之間的非線性耦合效應，還會與周圍結構發(fā)生流固耦合作用，其過程更加復雜，將在后續(xù)計劃中研究。