水下爆炸非均熵二維定常流的三族特征線解法*

李曉杰，楊晨琛，張程嬌，閆鴻浩，王小紅

(大連理工大學工程力學系工業裝備結構分析國家重點實驗室，遼寧 大連 116023)

水下爆炸是水下爆破與拆除、水中兵器毀傷研究、水面艦艇抗爆設計的基礎問題。最初采用半經驗公式與實驗相結合的方式開展研究[1]，主要的實驗研究方法有高頻傳感器法[2-3]和高速攝影法[4-5]。隨著計算機技術的飛速發展，數值模擬作為一種低成本且具有指導意義的研究方法，逐漸成為水下爆炸研究的重要輔助手段，其中商業有限元軟件在水下爆炸模擬中得到了廣泛應用[6-11]，一些新的算法如高精度差分方法[12-13]、無網格方法[14-15]、水平集方法[16-17]等不斷涌現。

水下爆炸沖擊波問題實際上是可壓縮流問題，其基本控制方程是雙曲型偏微分方程[18]，數值求解此類方程的離散方法包括有限差分法(finite difference method, FDM)、有限體積法、有限元法等。在沖擊波面上大多數離散方法采用人工黏性進行處理[19]，即在控制方程中添加一次、二次人工黏性項來光滑沖擊波面的強間斷，使間斷問題轉化為偏微分控制方程組能夠求解的連續問題，以引入人為誤差的代價抑制沖擊波面的數值跳躍。然而，人工黏性的使用在計算中存在問題：一是對于具體計算問題，選擇合適的人工黏性系數需要一定的經驗和反復嘗試；二是即便采用人工黏性也無法完全消除沖擊波后壓力等物理量的波動；三是對于固定的人工黏性系數取值，往往無法同時保證高、低壓段的精度，使得人為誤差在沖擊波較遠端和爆轟中心附近較為突出[8]。

特征線差分法[20]是一種將雙曲型偏微分方程組轉化為常微分方程組后計算流場的有限差分法，早期廣泛用于風洞、噴管等設計，在計算航天飛機外形的繞流問題上與實驗符合很好[21]。特征線在物理上是小擾動的傳播線，物理意義十分明確。特征線差分法實際上是求解流動的常微分方程組，在數學解法上有明確的理論依據，同時計算效率很高。在特征線差分法中，不再需要引入人工黏性處理沖擊波間斷，它將質點之間的相互作用看作沿特征線的擾動疊加，強間斷邊界只代表強一些的擾動，可以直接用強間斷方程確定空間位置以及處理壓力、流速、熵等物理量的跳躍。本文中從推導非等熵流的二維特征線方程出發，通過補充流線方程作為第三族特征線方程，得到適用求解非等熵流場的特征線差分方法，并將其應用于柱狀裝藥水下爆炸問題的計算。

1 二維定常非等熵流的特征線方程

二維定常可壓縮無黏流在歐拉坐標系下的控制方程如下[22]：

式中：ρ、p分別為流體密度和壓力；x、y為空間坐標；ux、uy為流體速度在坐標方向的投影；對于平面流動，δ取0，對于軸對稱流動，δ取1。考慮到壓力p和比內能e都是狀態量，可用另外兩個狀態量表示[23]，如寫成密度ρ和比熵s的函數p(ρ,s)和e(ρ,s)，得到dp/dρ的全導數形式：

式中：c為流體聲速，T為溫度。式(4)沿流體質點的跡線成立(在本文中與流線重合)，含義為壓力隨體導數與密度隨體導數的比例關系，而式(1)的最后一項代表密度隨體導數，聯立式(1)、式(4)并整理得到：

引入待定系數λ1、λ2，將式(2)、式(3)和式(5)線性相加，構造式(5)+式(2)·λ1ux+式(3)·λ2uy，得到新的方程：

假設存在方向導數：

式(7)實際上是含待定系數λ1、λ2的二元一次方程組，用流動偏轉角θ、流體馬赫角μ和流速u替換式(7)中u的x、y方向分量(ux=ucosθ,uy=usinθ)以及聲速(c=usinμ)，可得3組解(最后一組相當于平凡解)：

將這3組解代回式(7)和式(8)，最終得到3組方程：

即該極限趨近于一個有限值，以往做法是取軸線附近一點的值作為替代[21-22]，實際上取軸線上的du/2udx即可。

可以看出，式(10)～(12)所示的常微分方程組構成了二維定常可壓縮無黏非等熵流的完備三族特征線方程組，在形式上與等熵流特征線方程組[21]完全一致。其中非等熵的影響體現在相容方程包含的沿流線的熵變項中，相比于以往推導的非等熵特征線方程組[24]中所指的沿特征線的熵變項，顯然沿流線的熵變項才是非等熵流的物理實質。通過合理描述該熵變項，可求解帶化學反應或熱交換的非等熵流動，也可求解含沖擊波的非均熵流(沿流線等熵而流線間存在熵差)，但后者仍需第三族特征線(流線)確定溫度、比焓等熱力學量。只有對于沿流線沒有熵變且流線之間沒有熵差的均熵流，三族特征線才能退化為兩族特征線方程求解。

2 柱狀裝藥的水下爆炸模型

考慮柱狀裝藥水下爆炸過程，設裝藥無限長，爆轟穩定，爆速恒定不變，CJ(Chapman-Jouguet)面垂直于軸線的平面，忽略熱傳遞和界面失穩。若將參考系取在CJ面上，則可建立如圖1所示的定常模型[25]。爆轟產物的流動是二維定常可壓縮均熵無黏柱狀流，其中CJ面的馬赫數為1[26]。水中流動是二維定常可壓縮非均熵無黏柱狀流，水中沖擊波是來流速度為爆速的駐定曲沖擊波，沖擊波強度隨著徑向距離的增加而衰減，波前后產生熵跳躍，之后沿流線等熵，流線之間存在熵差，為非均熵流。在爆轟產物與水流的水-氣界面上，滿足壓力連續和法向速度連續的邊界條件，后者即為水、氣兩側偏轉角θ相等。在緊鄰爆轟波的水-氣界面處是曲沖擊波的起點，爆轟產物一側可以用普朗特-邁耶爾(Prandtl-Meyer)繞流描述[27]：

式中：Ma=u/c為爆轟產物的馬赫數。由于爆轟產物為均熵流，因此僅在流線上滿足的式(12)轉化為全流場適用的式(16)。再將水中斜沖擊波關系與爆轟產物普朗特-邁耶爾繞流方程聯立，即可求出流場的初始值。

具體計算時，對炸藥的爆轟產物選用JWL(Jones-Wilkins-Lee)狀態方程描述[28]：

式中：v0為炸藥初始比容，V=v/v0為相對體積；E為體積內能；A、B、R1、R2、ω為常數，且(?E/?p)ρ=V/ω。對水選用如下形式的Mie-Grüneison狀態方程描述：

式中：ρ0為水的初始密度，μ為壓縮率，μ=ρ/ρ0-1；e為比內能，A1、A2、A3、T1、T2、B0、B1均為系數，且有：

狀態方程的具體參數見表1和表2，其中炸藥PETN、TNT以及水的參數源于文獻[29]，SEP炸藥(其中PETN和石蠟的質量分數分別為65%和35%)的參數源于文獻[4]，D為爆速，pCJ為CJ壓力。在等熵卸載過程中，爆轟產物或水需要使用等熵條件：

表2 水的Mie-Grüneison狀態方程Table 2 Mie-Grüneison equation of state of water

3 水下爆炸計算實例

3.1 計算格式

非等熵流的三族特征線方程組實際上就是非齊次常微分方程組。以4節點計算格式(見圖2(a))為例，計算格式為：

式中：“[]Ⅰ,Ⅱ”代表沿特征線的加權取值；i、j為索引序號，代表邏輯空間的網格。為了避免隱式計算，常用顯格式的歐拉預估-校正法(經典特征線理論稱參數平均法[30])，即第一步預估的“[]Ⅰ,Ⅱ”取已知點的參數計算解點，第二步校正的“[]Ⅰ,Ⅱ”取已知點和解點參數的平均值，預估-校正法在形式上具有二階精度(Δx2,Δy2)。更普遍的5節點計算格式(見圖5(b))用4節點的計算格式通過內插獲得。

3.2 計算結果

在特征線差分法中，曲沖擊波的形狀通過一系列斜沖擊波逼近，逼近程度與差分步長相關。沖擊波形狀是壓力、速度等發生突躍變化的不連續面，用有限元軟件如AUTODYN也能大致觀測到。圖3展示了4種炸藥的特征線差分法計算結果(Cal., PETN)和AUTODYN模擬結果(計算區域1 500 mm×600 mm，網格尺寸1 mm×1 mm)，其中R、x為徑向和軸向距離，R0為裝藥半徑。圖4展示了SEP炸藥的特征線差分法計算結果和實驗結果[4]。可以看出：特征線差分法比較準確地捕捉到了近場沖擊波形狀。

根據斜沖擊波理論[21]，斜沖擊波上的物理參數如壓力、速度、偏轉角、激波角等，只有一個參數是獨立的。如果已知曲沖擊波形狀，那么從數學上可以確定曲沖擊波面上每一處的斜率(即激波角)，從而確定整個曲沖擊波面上的壓力、速度、偏轉角等。如果特征線差分法計算的近場沖擊波形狀足夠精確，可以推斷：特征線差分法計算的近場沖擊波壓力也足夠精確。圖5展示了特征線差分法和AUTODYN所代表的有限元法計算的某空間位置(R=2R0)的壓力時程曲線。可以看出：有限元法得到的壓力時程曲線具有明顯波動，而特征線差分法獲得的壓力時程曲線可以瞬間跳躍、連續衰減而無任何波動。

4 結 論

根據二維定常可壓縮超聲速非等熵流的控制方程，推導出其特征線方程，并通過補充流線方程作為第三族特征線方程，使特征線差分方法可以求解非等熵流問題，相比于以往添加沿特征線的熵變項，更能體現非等熵流的物理實質。

采用三族特征線方法，將非等熵流問題轉化為求解非齊次常微分方程組問題，提出了五點差分格式及其歐拉預估-校正解法，在理論上可以保證二階計算精度。

對柱狀裝藥的水下爆炸建立了定常模型，將非等熵特征線差分法應用于求解水下爆炸近場的非均熵流動。對幾種炸藥的水下爆炸近場沖擊波進行計算，結果表明，新的特征線差分方法可以準確地捕捉沖擊波形狀并獲得波后沖擊波壓力歷程，避免了常規方法計算沖擊波壓力的人工黏性誤差和數值波動，說明所提出的三族特征線差分法可以用于處理水下爆炸這種非均熵問題。

參考文獻:

[1] COLE R H. Underwater explosions[M]. New Jersy: Princeton University Press, 1948:3-13.

[2] STEINBERG D J. Spherical explosions and the equation of state of water: UCID-20974[R].Livermore, USA: Lawrence Livermore National Laboratory, 1987.

[3] 池家春,馬冰.TNT/RDX(40/60)炸藥球水中爆炸波研究[J].高壓物理學報,1999,13(3):199-204.

CHI Jiachun, MA Bing. Underwater explosion wave by a spherical charge of composition B-3[J].Chinese Journal of High Pressure Physics, 1999,13(3):199-204.

[4] ITOH S, SUZUKI O, NAGANO S, et al. Investigations of fundamental properties of underwater shock waves by high-speed photography[C]∥ 21st International Congress on: High-Speed Photography and Photonics. International Society for Optics and Photonics, 1995:916-928.

[5] 趙繼波,譚多望,李金河,等.TNT藥柱水中爆炸近場壓力軸向衰減規律[J].爆炸與沖擊,2008,28(6):539-543.

ZHAO Jibo, TAN Duowang, LI Jinhe, et al. Axial pressure damping of cylindrical TNT charges in the near underwater-explosion field[J]. Explosion and Shock Waves, 2008,28(6):539-543.

[6] 張振華,朱錫,白雪飛.水下爆炸沖擊波的數值模擬研究[J].爆炸與沖擊,2004,24(2):182-188.

ZHANG Zhenhua, ZHU Xi, BAI Xuefei. The study on numerical simulation of underwater blast wave[J]. Explosion and Shock Waves, 2004,24(2):182-188.

[7] 方斌,朱錫,張振華,等.水下爆炸沖擊波數值模擬中的參數影響[J].哈爾濱工程大學學報,2005,26(4):419-424.

FANG Bin, ZHU Xi, ZHANG Zhenhua, et al. Effect of parameters in numerical simulation of underwater shock wave[J]. Journal of Harbin Engineering University, 2005,26(4):419-424.

[8] 胡毅亭,賈憲振,饒國寧,等.水下爆炸沖擊波和氣泡脈動的數值模擬研究[J].艦船科學技術,2009,31(2):134-140.

HU Yiting, JIA Xianzhen, RAO Guoning, et al. Numerical study of underwater explosion shock wave and bubble pulse[J]. Ship Science and Technology, 2009,31(2):134-140.

[9] 劉科種,徐更光,辛春亮,等.AUTODYN水下爆炸數值模擬研究[J].爆破,2009,26(3):18-21.

LIU Kezhong, XU Gengguang, XIN Chunliang, et al. Research on numerical simulation in underwater explosion by AUTODYN[J]. Blasting, 2009,26(3):18-21.

[10] 辛春亮,秦健,劉科種,等.基于LS-DYNA軟件的水下爆炸數值模擬研究[J].彈箭與制導學報,2008,28(3):156-158.

XIN Chunliang, QIN Jian, LIU Kezhong, et al. Research on UNDEX numerical simulation based on LS-DYNA[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2008,28(3):156-158.

[11] 汪俊,劉建湖,李玉節.加筋圓柱殼水下爆炸動響應數值模擬[J].船舶力學,2006,10(2):126-137.

WANG Jun, LIU Jianhu, LI Yujie. Numerical simulation of dynamic response of ring-stiffened cylindrical shell subjected to underwater explosions[J]. Journal of Ship Mechanics, 2006,10(2):126-137.

[12] JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996,126(1):202-228.

[13] ZHANG X, SHU C W. On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[J]. Journal of Computational Physics, 2010,229(23):8918-8934.

[14] LIU M B, LIU G R, LAM K Y, et al. Smoothed particle hydrodynamics for numerical simulation of underwater explosion[J]. Computational Mechanics, 2003,30(2):106-118.

[15] ZHANG A, YANG W S, HUANG C, et al. Numerical simulation of column charge underwater explosion based on SPH and BEM combination[J]. Computers & Fluids, 2013,71:169-178.

[16] SUSSMAN M, SMEREKA P, OSHER S. A level set approach for computing solutions to incompressible two-phase flow[J]. Journal of Computational Physics, 1994,114(1):146-159.

[17] 師華強,汪玉,宗智,等.近水面水下爆炸二維Level-set數值模擬[J].計算力學學報,2010,27(3):409-414.

SHI Huaqiang, WANG Yu, ZONG Zhi, et al. 2D numerical simulation of underwater explosion near free surface based on level-set method[J]. Chinese Journal of Computational Mechanics, 2010,27(3):409-414.

[18] ANDERSON J D. Modern compressible flow with historical perspective[M]. 2nd ed. New York: McGraw-Hill, 1990:307-311.

[19] FLETCHER C A. Computational techniques for fluid dynamics[M]. Springer-Verlag, 1992:280-281.

[20] CHOU P C, HUANG S L, KARPP R R. Numerical calculation of blast waves by the method of characteristics[J]. AIAA Journal, 1967,5(4):618-623.

[21] RAKICH J, KUTLER P. Comparison of characteristics and shock capturing methods with application to the space shuttle vehicle[C]∥10th Aerospace Sciences Meeting, 1972:191.

[22] ZUCROW M J, HOFFMAN J D. 氣體動力學(上冊)[M].王汝涌,吳宗真,吳宗善,等譯.北京:國防工業出版社,1984:431-432.

[23] 李曉杰,張程嬌,閆鴻浩,等.水下爆炸近場非均熵流的特征線差分解法[J].爆炸與沖擊,2012,32(6):604-608.

LI Xiaojie, ZHANG Chengjiao, YAN Honghao, et al. Difference method of characteristics in isentropic flow of underwater explosion in near-field region[J]. Explosion and Shock Waves, 2012,32(6):604-608.

[24] HARTREE D R. Some practical methods of using characteristics in the calculation of non-steady compressible flow: AECU-2713[R]. Cambridge: Harvard University, 1953.

[25] STEBNOVSKII S V, CHERNOBAEV N N. Initial stage of an underwater explosion of cylindrical charges with foliated cases[J]. Combustion, Explosion and Shock Waves, 1982,18(3):358-362.

[26] 李曉杰,王金相,閆鴻浩,等.基于Laval噴管的圓管爆轟驅動近似解[J].工程爆破,2003,9(4):7-9.

LI Xiaojie, WANG Jinxiang, YAN Honghao, et al. The approximate solution of pipe driven by detonation based on the theory of de laval nozzle[J]. Engineering Blasting, 2003,9(4):7-9.

[27] 李曉杰,趙春風.基于通用炸藥狀態方程分析飛板運動規律的特征線法[J].爆炸與沖擊,2012,32(3):237-242.

LI Xiaojie, ZHAO Chunfeng. Characteristic curve method for analyzing movement of flyer plate based on universal equation of state of explosive[J]. Explosion and Shock Waves, 2012,32(3):237-242.

[28] 周清,張云海,毛允波.典型民用建筑內爆炸波傳播機理[J].山東科技大學學報(自然科學版),2009,28(6):36-44.

ZHOU Qing, ZHANG Yunhai, MAO Yunbo. The spreading mechanism of explosion waves in typical civil building[J]. Journal of Shandong University of Science and Technology (Natural Science), 2009,28(6):36-44.

[29] ANSYS Inc. ANSYS Autodyn user manual, release 13.0[M].ANSYS Inc., 2010.

[30] ZUCROW M J, HOFFMAN J D. 氣體動力學(下冊)[M].魏叔如,吳宗善,王汝涌,等譯.北京:國防工業出版社,1984:143-144.