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

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

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

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

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

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

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

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

式中：ρ、p分別為流體密度和壓力；x、y為空間坐標(biāo)；ux、uy為流體速度在坐標(biāo)方向的投影；對(duì)于平面流動(dòng)，δ取0，對(duì)于軸對(duì)稱流動(dòng)，δ取1?？紤]到壓力p和比內(nèi)能e都是狀態(tài)量，可用另外兩個(gè)狀態(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的二元一次方程組，用流動(dòng)偏轉(zhuǎn)角θ、流體馬赫角μ和流速u替換式(7)中u的x、y方向分量(ux=ucosθ,uy=usinθ)以及聲速(c=usinμ)，可得3組解(最后一組相當(dāng)于平凡解)：

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

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

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

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

考慮柱狀裝藥水下爆炸過(guò)程，設(shè)裝藥無(wú)限長(zhǎng)，爆轟穩(wěn)定，爆速恒定不變，CJ(Chapman-Jouguet)面垂直于軸線的平面，忽略熱傳遞和界面失穩(wěn)。若將參考系取在CJ面上，則可建立如圖1所示的定常模型[25]。爆轟產(chǎn)物的流動(dòng)是二維定?？蓧嚎s均熵?zé)o黏柱狀流，其中CJ面的馬赫數(shù)為1[26]。水中流動(dòng)是二維定?？蓧嚎s非均熵?zé)o黏柱狀流，水中沖擊波是來(lái)流速度為爆速的駐定曲沖擊波，沖擊波強(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)化為全流場(chǎng)適用的式(16)。再將水中斜沖擊波關(guān)系與爆轟產(chǎn)物普朗特-邁耶爾繞流方程聯(lián)立，即可求出流場(chǎng)的初始值。

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

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

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

狀態(tài)方程的具體參數(shù)見(jiàn)表1和表2，其中炸藥PETN、TNT以及水的參數(shù)源于文獻(xiàn)[29]，SEP炸藥(其中PETN和石蠟的質(zhì)量分?jǐn)?shù)分別為65%和35%)的參數(shù)源于文獻(xiàn)[4]，D為爆速，pCJ為CJ壓力。在等熵卸載過(guò)程中，爆轟產(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ì)算格式(見(jiàn)圖2(a))為例，計(jì)算格式為：

式中：“[]Ⅰ,Ⅱ”代表沿特征線的加權(quán)取值；i、j為索引序號(hào)，代表邏輯空間的網(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ì)算格式(見(jiàn)圖5(b))用4節(jié)點(diǎn)的計(jì)算格式通過(guò)內(nèi)插獲得。

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

在特征線差分法中，曲沖擊波的形狀通過(guò)一系列斜沖擊波逼近，逼近程度與差分步長(zhǎng)相關(guān)。沖擊波形狀是壓力、速度等發(fā)生突躍變化的不連續(xù)面，用有限元軟件如AUTODYN也能大致觀測(cè)到。圖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)確地捕捉到了近場(chǎng)沖擊波形狀。

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

4 結(jié) 論

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

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

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

參考文獻(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藥柱水中爆炸近場(chǎng)壓力軸向衰減規(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] 胡毅亭,賈憲振,饒國(guó)寧,等.水下爆炸沖擊波和氣泡脈動(dòng)的數(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é).加筋圓柱殼水下爆炸動(dòng)響應(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. 氣體動(dòng)力學(xué)(上冊(cè))[M].王汝涌,吳宗真,吳宗善,等譯.北京:國(guó)防工業(yè)出版社,1984:431-432.

[23] 李曉杰,張程嬌,閆鴻浩,等.水下爆炸近場(chǎng)非均熵流的特征線差分解法[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ū)動(dòng)近似解[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)動(dòng)規(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. 氣體動(dòng)力學(xué)(下冊(cè))[M].魏叔如,吳宗善,王汝涌,等譯.北京:國(guó)防工業(yè)出版社,1984:143-144.