李曉杰，張程嬌，閆鴻浩，王小紅，王宇新

（大連理工大學工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連116024）

水下爆炸的過程分為炸藥的爆轟、沖擊波在水中的形成及傳播、氣體球膨脹收縮運動等3個部分［1］。以往對遠場中低壓力區(qū)的爆炸沖擊波和氣泡脈動的研究較多，其中遠場沖擊波壓力與氣泡脈動壓力的測試常被用于炸藥爆炸能量評定。對水下爆炸沖擊波傳播問題的研究通常采用實驗和數(shù)值模擬的方法。實驗測試近場沖擊波時，由于沖擊波峰值壓力可以達到數(shù)萬至數(shù)十萬個大氣壓，沖擊波對壓力傳感器性能要求很高并容易造成嚴重破壞使得沖擊波壓力難于測量，因此水下爆炸近場沖擊波的實驗研究比較困難。近年來對水下爆炸近場問題的研究已逐步采用數(shù)值計算的方法解決。目前采用的數(shù)值計算方法很多，但本質(zhì)上仍然是有限差分和有限元2類，采用有限元或高精差分方法模擬沖擊波間斷問題時，在計算中必須引入人工粘性，對沖擊波間斷做連續(xù)化處理，所計算的沖擊波壓力受人工粘性系數(shù)取值的影響［2－4］。盡管人們已經(jīng)通過將大量計算結(jié)果與實驗比較，獲得了較合適的人工粘性系數(shù)，但人工粘性系數(shù)的使用引入了人為干預因素。

特征線差分方法［5］是早期流場計算中常用的一種方法，特征線差分法具有物理概念明確、計算簡單、計算精度高等特點。但由于特征線差分法是根據(jù)流場全域等熵假定（即均熵假定）獲得的，所以只能用于研究均熵流動問題。盡管早期學者也使用特征線法對水下爆炸問題進行過研究，但一直忽略了水下變強度沖擊波后的非均熵流特征，而是采用簡單的均熵假設來近似處理沖擊波問題，這使得特征線差分法計算的正確性受到了質(zhì)疑。但是特征線差分方法的優(yōu)點非常明顯：（1）無需對沖擊波做連續(xù)化處理，（2）可以直接將流體的偏微分方程轉(zhuǎn)化為計算精度較高的常微分方程來求解，（3）特征線本身對應著聲波的傳輸軌跡具有明確的物理意義；所以該方法的研究一直受到重視。

本文中從可壓縮流的基本方程出發(fā)，重新推導含有熵變修正項的2族特征線方程，然后將質(zhì)點的跡線方程及其能量方程補充為第3族特征線方程，使得特征線方法適用于非均熵流場的求解。而后利用所獲得的3族特征線方程組，對質(zhì)量分數(shù)比w（TNT）∶w（RDX）＝40∶60的TNT／RDX球形炸藥水下爆炸進行數(shù)值計算，并將沖擊波近場壓力的計算結(jié)果與文獻［6］中的實驗數(shù)據(jù)進行比較，驗證這種改進的特征線差分法的正確性。

1 非均熵流的特征線差分方程

1.1 非均熵流的特征線方程推導

水下爆炸作用時間極短，通常忽略傳熱效應，視為絕熱過程。忽略熱傳導效應和粘性效應的一維可壓縮流的連續(xù)方程和運動方程如下

式中：ρ、u、p、r和t分別為流體密度、速度、壓強、空間坐標和時間坐標。N＝0，1，2分別表示一維平面、柱面和球面流動問題。以上2式是拉格朗日形式的連續(xù)方程和運動方程，也就是沿質(zhì)點跡線的流動方程。流場內(nèi)質(zhì)點的壓力可以用熵S與密度ρ來表示p＝p（ρ，S），可以寫成

式中：下標III表示物理量沿質(zhì)點跡線的變化關系，引入系數(shù)λ，聯(lián)立式（4）與式（2）得到

若令λ＝±ρc，那么上式就可以寫成

若時空平面內(nèi)的2族特征線方程為dr／dt＝u±c，代入式（6）就可以獲得物理平面上的相容關系，其中右行特征線方程與相容關系如下

左行特征線方程與相容關系如下

式中：（dQ）III＝（TdS）III，能量方程中顯含有熵變項，對于一般問題必須給出熵S的具體表達式。但對于水下爆炸這類高速沖擊問題，依然可以假設流體質(zhì)點之間是絕熱的，即每個質(zhì)點沿跡線是等熵的，但各質(zhì)點之間的熵變互不相同。因此不要求全流場具有相同的熵變，這就是熵變沿流體質(zhì)點跡線等熵的全場非均熵流，即（dS）III＝0，因此能量方程就可以寫成隱含熵S的形式，即

同樣地，由于（dS）III＝0，式（7）、（8）中的項也為0。式（7）、（8）、（10）就構(gòu)成了可以用于求解非均熵可壓縮流的隱含熵S項3族特征線方程。應該注意的是，盡管這3族特征線方程與以往均熵流的方程在形式上幾乎完全相同，但所代表的物理意義卻完全不同，可以用于處理初始熵變非均勻分布的絕熱非均熵流了。例如：對于球?qū)ΨQ一維水下爆炸這類變沖擊波強度的問題，由沖擊波強度確定了波陣面的熵變后，沖擊波后方的流場就變成這種沿質(zhì)點跡線等熵的全場非均熵流場；另外，對于單個沖擊波追趕卸載問題也具有同樣的特征，也可以化為這類問題用3族特征線方程組求解。

從上述方程可見，在時空平面內(nèi)，第I、II族特征線方程與標準特征線方程完全相同，均為dr／dt＝u±c；在物理平面內(nèi)，第I、II族特征線代表了聲波小擾動的傳播線，推導的特征線方程和標準特征線方程具有相同的物理意義，不同在于推導的相容關系中包含了熵變影響項。這也說明了在任何復雜的流場中，其中任意一點的擾動均會沿著聲波小擾動線傳播，即第I、II族特征線是物理存在的，而復雜問題在于如何確定特征線上的相容關系。第3族特征線方程組（10）采用了絕熱等熵假設，如果使用更復雜的能量方程，該方法可以用于求解帶化學反應或帶耗散的復雜流動。

1.2 非均熵流的特征線差分格式

在應用標準特征線法求解均熵流場的問題時，只需要前2族特征線就可以求解，其求解方法是用一條右行特征線I與一條左行特征線II相交于某一點，該點的各個參數(shù)由這2條特征線所通過的上一時刻的2個節(jié)點參數(shù)來確定。如圖1所示，右行特征線I與左行特征線II相交于A點，A點的各個參數(shù)可以用這2條特征線所通過的前一時刻的B、C等2點來求解；如此在給定了計算域初值條件的情況下，就可以求解出域內(nèi)相關點的參數(shù)。標準特征線法求解過程簡單并且物理意義明確，然而這種方法是以計算域內(nèi)流場均熵為基礎的，即對均熵流場求解，在計算過程中并未考慮相鄰質(zhì)點間熵值不相等對計算結(jié)果的影響。

圖1 標準特征線法在時空域內(nèi)的差分方法示意圖Fig.1 Standard difference method of characteristics in r－t plane

圖2 改進的特征線法的差分格式Fig.2 Difference scheme for improved difference method of characteristics

在絕熱的流場中，各質(zhì)點的熵值沿跡線保持不變，本文的特征線方程中添加了熵修正項，并補充加入了質(zhì)點的跡線方程作為第3族特征線方程，因此可以求解非均熵流場。由于引入了3族特征線也就再無法使用簡單的差分方法進行計算了，因此必須使用新的差分方法進行求解。改進特征線法的差分求解方法如下：如圖2所示，假定在t0時刻點A、B0、C的參數(shù)已知，經(jīng)過時間Δt，質(zhì)點由B0沿跡線III運動至B1處，在時空平面上可以沿跡線III求出B1點的坐標；通過B1點的2條特征線I和II又與t0時間線相交于點E和點F，這2點的坐標及各參數(shù)可以通過在A、B0與B0、C之間分別插值求得；最后再用點E和點F就可以求出B1點處質(zhì)點的各個狀態(tài)參數(shù)。若給定初值條件和邊界條件，則如此逐層迭代計算就可得出整個域內(nèi)相關點的參數(shù)。

2 水下爆炸近場壓力分布的數(shù)值計算

一維球面水下爆炸問題是典型的非均熵流問題。爆轟波入射到水中產(chǎn)生球面的水中沖擊波，而水中沖擊波強度隨傳播距離增加而不斷衰減。由于沖擊波引起的熵增與沖擊波強度有關，因此水下爆炸沖擊波后的熵增也隨著傳播距離增加而不斷減小，所以說球形炸藥水下爆炸是典型的非均熵流問題。

2.1 計算模型

將改進的特征線法程序化，對質(zhì)量分數(shù)比w（TNT）∶w（RDX）＝40∶60的TNT／RDX炸藥的水下爆炸問題進行數(shù)值計算。選用JWL狀態(tài)方程

式中：e為單位質(zhì)量內(nèi)能；V 為相對體積，V＝v／v0，v0＝1／ρ0，ρ0為炸藥初始密度，ρ0＝1.717g／cm3；方程中的其他參數(shù)分別取值為：A＝524.23GPa，B＝7.678GPa，R1＝4.2，R2＝1.1，ω＝0.34；炸藥爆速D＝7.98km／s。

對水選用多項式形式的狀態(tài)方程，當水受壓和膨脹時，有

式中：μ 是壓縮度，μ＝ρ／ρ0－1；e是單位質(zhì)量內(nèi)能；ρ0是初始密度，ρ0＝1.0g／cm3；方程中的其他參數(shù)：A1＝2.2GPa，A2＝9.54GPa，A3＝14.57GPa，T1＝2.2GPa，T2＝0，B0＝0.28，B1＝0.28。

計算中采用一維球?qū)ΨQ模型，在無限水域中放入1kg的質(zhì)量分數(shù)比w（TNT）∶w（RDX）＝40∶60的TNT／RDX炸藥球，藥球半徑為5.18cm。為了避開藥球中心處的計算奇點，起爆點設置在距離球心0.2cm的球面上，起爆后在該處設置固壁邊界條件；網(wǎng)格大小設置為0.249cm，炸藥內(nèi)劃分20個網(wǎng)格。

2.2 計算結(jié)果

數(shù)值計算結(jié)果如圖3所示，圖中繪制出了1≤R／R0≤4范圍內(nèi)的峰值壓力pmax－比例距離R／R0的曲線。為檢驗本文中提出的改進的特征線法，將計算結(jié)果與池家春等［6］的實驗結(jié)果進行了比較，實驗值公式如下

從圖3中可以看出，沖擊波峰值壓力的計算值與實驗值符合較好。在比例距離1≤R／R0≤4的近場范圍內(nèi)，兩者最大相對誤差為13.1%，這說明本文的特征線法完全可以用于處理水下爆炸近場的非均熵流問題。

圖3 近場沖擊波峰值壓力的比較Fig.3 Comparison of the peak pressures in the near－field region

3 結(jié) 論

本文中提出了一種改進的特征線法，在特征線方程中添加了熵的修正項，并將質(zhì)點的跡線作為第3族特征線加入到特征線方程組中，使它可以用于求解變強度沖擊波及波后的非均熵流場。由于時空平面內(nèi)的第I、II族特征線與標準特征線方程在形式上完全相同，因此非均熵流的前2族特征線方程同樣也代表了聲波小擾動的傳播線，與標準特征線方程具有相同的物理意義。

將改進的特征線法程序化，對質(zhì)量分數(shù)比w（TNT）∶w（RDX）＝40∶60的TNT／RDX炸藥球的水下爆炸問題做了計算。在比例距離1≤R／R0≤4的近場范圍內(nèi)，用本文的特征線法計算得到的沖擊波峰值壓力與實驗值［6］最大誤差為13.1%，從而驗證了非均熵特征線法的正確性。另外，這種方法無需采用人工粘性系數(shù)處理沖擊波間斷，只需給定計算所需的相關物理參數(shù)即可求解，因此減少了對沖擊波的人為干預因素。

［1］Cole R H.Underwater explosions［M］.New Jersy：Princeton University Press，1948：3－13，110.

［2］劉科種，徐更光，辛春亮，等.AUTODYN水下爆炸數(shù)值模擬研究［J］.爆破，2009，26（3）：18－21.LIU Ke－zhong，XU Geng－guang，XIN Chun－liang，et al.Research on numerical simulation in underwater explosioln by AUTODYN［J］.Blasting，2009，26（3）：18－21.

［3］徐豫新，王樹山，李園.水下爆炸數(shù)值仿真研究［J］.彈箭與制導學報，2009，29（6）：95－97，102.XU Yu－xin，WANG Shu－shan，LI Yuan.Study on numerical simulation of the underwater explosive［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2009，29（6）：95－97，102.

［4］方斌，朱錫，張振華，等.水下爆炸沖擊波數(shù)值模擬中的參數(shù)影響［J］.哈爾濱工程大學學報，2005，26（4）：419－424.FANG Bin，ZHU Xi，ZHANG Zhen－h(huán)ua，et al.Effect of parameters in numerical simulation of underwater shock wave［J］.Journal of Harbin Engineering University，2005，26（4）：419－424.

［5］Chou P C，Huang S L，Karpp R R.Numerical calculation of blast wave by the method of characteristics［J］.AIAA Journal，1967，5（4）：618－623.

［6］池家春，馬冰.TNT／RDX（40／60）炸藥球水中爆炸波研究［J］.高壓物理學報，1999，13（3）：199－204.CHI Jia－chun，MA Bing.Underwater explosion wave by a spherical charge of composition B－3［J］.Chinese Journal of High Pressure Physics，1999，13（3）：199－244.