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

李晓杰，杨晨琛，张程娇，闫鸿浩，王小红

(大连理工大学工程力学系工业装备结构分析国家重点实验室，辽宁 大连 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.