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

(大連理工大學(xué)工程力學(xué)系工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023)

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

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

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

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

二維定?？蓧嚎s無黏流在歐拉坐標(biāo)系下的控制方程如下[22]：

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

式中：c為流體聲速，T為溫度。式(4)沿流體質(zhì)點(diǎn)的跡線成立(在本文中與流線重合)，含義為壓力隨體導(dǎo)數(shù)與密度隨體導(dǎo)數(shù)的比例關(guān)系，而式(1)的最后一項(xiàng)代表密度隨體導(dǎo)數(shù)，聯(lián)立式(1)、式(4)并整理得到：

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

假設(shè)存在方向?qū)?shù)：

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

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

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

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

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

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

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

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

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

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

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

表2 水的Mie-Grüneison狀態(tài)方程Table 2 Mie-Grüneison equation of state of water

3 水下爆炸計(jì)算實(shí)例

3.1 計(jì)算格式

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

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

3.2 計(jì)算結(jié)果

在特征線差分法中，曲沖擊波的形狀通過一系列斜沖擊波逼近，逼近程度與差分步長相關(guān)。沖擊波形狀是壓力、速度等發(fā)生突躍變化的不連續(xù)面，用有限元軟件如AUTODYN也能大致觀測到。圖3展示了4種炸藥的特征線差分法計(jì)算結(jié)果(Cal., PETN)和AUTODYN模擬結(jié)果(計(jì)算區(qū)域1 500 mm×600 mm，網(wǎng)格尺寸1 mm×1 mm)，其中R、x為徑向和軸向距離，R0為裝藥半徑。圖4展示了SEP炸藥的特征線差分法計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果[4]?？梢钥闯觯禾卣骶€差分法比較準(zhǔn)確地捕捉到了近場沖擊波形狀。

根據(jù)斜沖擊波理論[21]，斜沖擊波上的物理參數(shù)如壓力、速度、偏轉(zhuǎn)角、激波角等，只有一個參數(shù)是獨(dú)立的。如果已知曲沖擊波形狀，那么從數(shù)學(xué)上可以確定曲沖擊波面上每一處的斜率(即激波角)，從而確定整個曲沖擊波面上的壓力、速度、偏轉(zhuǎn)角等。如果特征線差分法計(jì)算的近場沖擊波形狀足夠精確，可以推斷：特征線差分法計(jì)算的近場沖擊波壓力也足夠精確。圖5展示了特征線差分法和AUTODYN所代表的有限元法計(jì)算的某空間位置(R=2R0)的壓力時程曲線?？梢钥闯觯河邢拊ǖ玫降膲毫r程曲線具有明顯波動，而特征線差分法獲得的壓力時程曲線可以瞬間跳躍、連續(xù)衰減而無任何波動。

4 結(jié) 論

根據(jù)二維定常可壓縮超聲速非等熵流的控制方程，推導(dǎo)出其特征線方程，并通過補(bǔ)充流線方程作為第三族特征線方程，使特征線差分方法可以求解非等熵流問題，相比于以往添加沿特征線的熵變項(xiàng)，更能體現(xiàn)非等熵流的物理實(shí)質(zhì)。

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

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

參考文獻(xiàn):

[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].高壓物理學(xué)報(bào),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藥柱水中爆炸近場壓力軸向衰減規(guī)律[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] 張振華,朱錫,白雪飛.水下爆炸沖擊波的數(shù)值模擬研究[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] 方斌,朱錫,張振華,等.水下爆炸沖擊波數(shù)值模擬中的參數(shù)影響[J].哈爾濱工程大學(xué)學(xué)報(bào),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] 胡毅亭,賈憲振,饒國寧,等.水下爆炸沖擊波和氣泡脈動的數(shù)值模擬研究[J].艦船科學(xué)技術(shù),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水下爆炸數(shù)值模擬研究[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軟件的水下爆炸數(shù)值模擬研究[J].彈箭與制導(dǎo)學(xué)報(bào),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] 汪俊,劉建湖,李玉節(jié).加筋圓柱殼水下爆炸動響應(yīng)數(shù)值模擬[J].船舶力學(xué),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] 師華強(qiáng),汪玉,宗智,等.近水面水下爆炸二維Level-set數(shù)值模擬[J].計(jì)算力學(xué)學(xué)報(bào),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. 氣體動力學(xué)(上冊)[M].王汝涌,吳宗真,吳宗善,等譯.北京:國防工業(yè)出版社,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噴管的圓管爆轟驅(qū)動近似解[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] 李曉杰,趙春風(fēng).基于通用炸藥狀態(tài)方程分析飛板運(yùn)動規(guī)律的特征線法[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] 周清,張?jiān)坪?毛允波.典型民用建筑內(nèi)爆炸波傳播機(jī)理[J].山東科技大學(xué)學(xué)報(bào)(自然科學(xué)版),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. 氣體動力學(xué)(下冊)[M].魏叔如,吳宗善,王汝涌,等譯.北京:國防工業(yè)出版社,1984:143-144.