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

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

doi:10.11883/bzycj-2016-0314

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

doi: 10.11883/bzycj-2016-0314

大连理工大学工程力学系工业装备结构分析国家重点实验室, 辽宁大连 116023

基金项目:

国家自然科学基金项目 11272081

国家自然科学基金项目 11672067

详细信息

作者简介:
李晓杰(1963-), 男, 教授, 博士生导师, dalian03@qq.com

中图分类号: O381
计量
- 文章访问数: 5739
- HTML全文浏览量: 2168
- PDF下载量: 245
- 被引次数: 3
出版历程
- 收稿日期: 2016-10-14
- 修回日期: 2017-04-05
- 刊出日期: 2018-07-25

A finite difference method of three characteristic lines of two-dimensional non-isentropic steady flow of cylindrical explosive underwater explosion

State Key Laboratory of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics, Dalian University of Technology, Dalian 116023, Liaoning, China

摘要

摘要: 针对二维定常可压缩超声速非等熵柱状流，提出一种特征线差分解法，通过在沿马赫线的相容方程中添加沿流线的熵变项以描述非等熵效应，得到等熵流和非等熵流均适用的三族特征线方程组。根据水下爆炸近场特点，建立无限长柱状装药的定常模型，将三族特征线方程组用有限差分法离散求解，通过构造合适的网格保证计算格式可以数值上收敛，由此编制程序并计算几种柱状炸药的水下爆炸近场冲击波。对比有限元模拟结果和实验结果发现，特征线差分法可以比较准确地捕捉冲击波形状并计算冲击波后流场，从而验证了所提出的三族特征线差分法的准确性。
- 水下爆炸 /
- 柱状装药 /
- 非均熵流 /
- 近场冲击波 /
- 特征线差分法
Abstract: In the present paper, we proposed a two-dimensional finite difference method (FDM) of characteristic lines to address problems of the non-isentropic steady flow of cylindrical explosive underwater explosion. This method describes the non-isentropic effect by adding an entropy-related variable along the flow line to the pressure-related equation along the Mach line, so that both the isentropic flow and the non-isentropic flow can be described in the same equations of the characteristics. Based on the features of the near-field shock wave we firstly modeled the underwater explosion with an infinitely long cylindrical explosive, then discretized those equations using this finite difference method and constructed an appropriate grid to ensure the numerical convergence, and finally calculated the underwater near-field shock wave for several explosives by programming. The numerical examples showed that the results of this method are consistent with those of the commercial finite element software AUTODYN and those of experiments, suggesting that the FDM of characteristics can capture the shock wave front with relatively high accuracy, and confirming that this method is applicable to solving problems in cylindrical explosive underwater explosion.
- underwater explosion /
- cylindrical charge /
- non-isentropic flow /
- near-field shock wave /
- finite difference method of characteristics

HTML全文

水下爆炸是水下爆破与拆除、水中兵器毁伤研究、水面舰艇抗爆设计的基础问题。最初采用半经验公式与实验相结合的方式开展研究^[1]，主要的实验研究方法有高频传感器法^[2-3]和高速摄影法^[4-5]。随着计算机技术的飞速发展，数值模拟作为一种低成本且具有指导意义的研究方法，逐渐成为水下爆炸研究的重要辅助手段，其中商业有限元软件在水下爆炸模拟中得到了广泛应用^[6-11]，一些新的算法如高精度差分方法^[12-13]、无网格方法^[14-15]、水平集方法^[16-17]等不断涌现。

水下爆炸冲击波问题实际上是可压缩流问题，其基本控制方程是双曲型偏微分方程^[18]，数值求解此类方程的离散方法包括有限差分法(finite difference method, FDM)、有限体积法、有限元法等。在冲击波面上大多数离散方法采用人工黏性进行处理^[19]，即在控制方程中添加一次、二次人工黏性项来光滑冲击波面的强间断，使间断问题转化为偏微分控制方程组能够求解的连续问题，以引入人为误差的代价抑制冲击波面的数值跳跃。然而，人工黏性的使用在计算中存在问题：一是对于具体计算问题，选择合适的人工黏性系数需要一定的经验和反复尝试；二是即便采用人工黏性也无法完全消除冲击波后压力等物理量的波动；三是对于固定的人工黏性系数取值，往往无法同时保证高、低压段的精度，使得人为误差在冲击波较远端和爆轰中心附近较为突出^[8]。

特征线差分法^[20]是一种将双曲型偏微分方程组转化为常微分方程组后计算流场的有限差分法，早期广泛用于风洞、喷管等设计，在计算航天飞机外形的绕流问题上与实验符合很好^[21]。特征线在物理上是小扰动的传播线，物理意义十分明确。特征线差分法实际上是求解流动的常微分方程组，在数学解法上有明确的理论依据，同时计算效率很高。在特征线差分法中，不再需要引入人工黏性处理冲击波间断，它将质点之间的相互作用看作沿特征线的扰动叠加，强间断边界只代表强一些的扰动，可以直接用强间断方程确定空间位置以及处理压力、流速、熵等物理量的跳跃。本文中从推导非等熵流的二维特征线方程出发，通过补充流线方程作为第三族特征线方程，得到适用求解非等熵流场的特征线差分方法，并将其应用于柱状装药水下爆炸问题的计算。

1. 二维定常非等熵流的特征线方程

二维定常可压缩无黏流在欧拉坐标系下的控制方程如下^[22]：

$\frac{{\partial {u_x}}}{{\partial x}} + \frac{{\partial {u_y}}}{{\partial y}} + \frac{{\delta {u_y}}}{y} + \frac{1}{\rho }\left( {{u_x}\frac{{\partial \rho }}{{\partial x}} + {u_y}\frac{{\partial \rho }}{{\partial y}}} \right) = 0$

(1)

${u_x}\frac{{\partial {u_x}}}{{\partial x}} + {u_y}\frac{{\partial {u_x}}}{{\partial y}} + \frac{1}{\rho }\frac{{\partial \rho }}{{\partial x}} = 0$

(2)

${u_x}\frac{{\partial {u_y}}}{{\partial x}} + {u_y}\frac{{\partial {u_y}}}{{\partial y}} + \frac{1}{\rho }\frac{{\partial \rho }}{{\partial y}} = 0$

(3)

式中：ρ、p分别为流体密度和压力；x、y为空间坐标；u_x、u_y为流体速度在坐标方向的投影；对于平面流动，δ取0，对于轴对称流动，δ取1。考虑到压力p和比内能e都是状态量，可用另外两个状态量表示^[23]，如写成密度ρ和比熵s的函数p(ρ, s)和e(ρ, s)，得到dp/dρ的全导数形式：

$\frac{{{\rm{d}}p}}{{{\rm{d}}\rho }} = {\left( {\frac{{\partial p}}{{\partial \rho }}} \right)_s} + {\left( {\frac{{\partial p}}{{\partial s}}} \right)_\rho }\frac{{{\rm{d}}s}}{{{\rm{d}}\rho }} = {c^2} + {\left( {\frac{{\partial p}}{{\partial e}}} \right)_\rho }{\left( {\frac{{\partial e}}{{\partial s}}} \right)_\rho }\frac{{{\rm{d}}s}}{{{\rm{d}}\rho }} = {c^2} + {\left( {\frac{{\partial p}}{{\partial e}}} \right)_\rho }\frac{{T{\rm{d}}s}}{{{\rm{d}}\rho }}$

(4)

式中：c为流体声速，T为温度。式(4)沿流体质点的迹线成立(在本文中与流线重合)，含义为压力随体导数与密度随体导数的比例关系，而式(1)的最后一项代表密度随体导数，联立式(1)、式(4)并整理得到：

${c^2}\frac{{\partial {u_x}}}{{\partial x}} + {c^2}\frac{{\partial {u_y}}}{{\partial y}} + \frac{1}{\rho }\left( {{u_x}\frac{{\partial p}}{{\partial x}} + {u_y}\frac{{\partial p}}{{\partial y}}} \right) - \\\frac{1}{\rho }\left( {{u_x}\frac{{\partial \rho }}{{\partial x}} + {u_y}\frac{{\partial \rho }}{{\partial y}}} \right){\left( {\frac{{\partial p}}{{\partial e}}} \right)_\rho }\frac{{T{\rm{d}}s}}{{{\rm{d}}\rho }} + \frac{{\delta {c^2}{u_y}}}{y} = 0$

(5)

引入待定系数λ₁、λ₂，将式(2)、式(3)和式(5)线性相加，构造式(5)+式(2)·λ₁u_x+式(3)·λ₂u_y，得到新的方程：

$\begin{array}{l} \left( {{c^2} + {\lambda _1}u_x^2} \right)\left( {\frac{{\partial {u_x}}}{{\partial x}} + \frac{{{\lambda _1}{u_x}{u_y}}}{{{c^2} + {\lambda _1}u_x^2}}\frac{{\partial {u_x}}}{{\partial y}}} \right) + {\lambda _2}{u_x}{u_y}\left( {\frac{{\partial {u_y}}}{{\partial x}} + \frac{{{c^2} + {\lambda _2}u_y^2}}{{{\lambda _2}{u_x}{u_y}}}\frac{{\partial {u_y}}}{{\partial y}}} \right) + \\ \frac{{\left( {1 + {\lambda _1}} \right){u_x}}}{\rho }\left[ {\frac{{\partial p}}{{\partial x}} + \frac{{\left( {1 + {\lambda _2}} \right){u_y}}}{{\left( {1 + {\lambda _1}} \right){u_x}}}\frac{{\partial p}}{{\partial y}}} \right] - \frac{{{u_x}}}{\rho }\left( {\frac{{\partial \rho }}{{\partial x}} + \frac{{{u_y}}}{{{u_x}}}\frac{{\partial \rho }}{{\partial y}}} \right){\left( {\frac{{\partial p}}{{\partial e}}} \right)_\rho }\frac{{T{\rm{d}}s}}{{{\rm{d}}\rho }} + \frac{{\delta {c^2}{u_y}}}{y} = 0 \end{array}$

(6)

假设存在方向导数：

$\frac{{{\rm{d}}y}}{{{\rm{d}}x}}\left| {_{{\rm{dir}}}} \right. = \frac{{{\lambda _1}{u_x}{u_y}}}{{{c^2} + {\lambda _1}u_x^2}} \equiv \frac{{{c^2} + {\lambda _2}u_y^2}}{{{\lambda _2}{u_x}{u_y}}} \equiv \frac{{\left( {1 + {\lambda _2}} \right){u_y}}}{{\left( {1 + {\lambda _1}} \right){u_x}}}$

(7)

将式(7)代入式(6)，考虑到 $\frac{\partial \rho }{\partial x}+\frac{{{u}_{y}}}{{{u}_{x}}}\frac{\partial \rho }{\partial y}$ 表示密度ρ沿流线的导数 $\frac{\text{d}\rho }{\text{d}x}\left| _{\text{stream}} \right.$ ，因此 $\left( \frac{\partial \rho }{\partial x}\text{+}\frac{{{u}_{y}}}{{{u}_{x}}}\frac{\partial \rho }{\partial y} \right)\frac{{{\text{d}}_{s}}}{\text{d}\rho }$ 也就表示比熵s沿流线的导数 $\frac{\text{d}s}{\text{d}x}\left| _{\text{stream}} \right.$ ，最终式(6)变换为如下方程：

$\left( {{c^2} + {\lambda _1}u_x^2} \right)\frac{{{\rm{d}}{u_x}}}{{{\rm{d}}x}}\left| {_{{\rm{dir}}}} \right. + {\lambda _2}{u_x}{u_y}\frac{{{\rm{d}}{u_y}}}{{{\rm{d}}x}}\left| {_{{\rm{dir}}}} \right. + \\\frac{{\left( {1 + {\lambda _1}} \right){u_x}}}{\rho }\frac{{{\rm{d}}p}}{{{\rm{d}}x}}\left| {_{{\rm{dir}}}} \right. - \frac{{{u_x}}}{\rho }{\left( {\frac{{\partial p}}{{\partial e}}} \right)_\rho }\frac{{T{\rm{d}}s}}{{{\rm{d}}x}}\left| {_{{\rm{stream}}}} \right. + \frac{{\delta {c^2}{u_y}}}{y} = 0$

(8)

式(7)实际上是含待定系数λ₁、λ₂的二元一次方程组，用流动偏转角θ、流体马赫角μ和流速u替换式(7)中u的x、y方向分量(u_x=u cos θ, u_y=u sin θ)以及声速(c=u sin μ)，可得3组解(最后一组相当于平凡解)：

$\begin{array}{l} {\rm{I}}:\;\;\;{\lambda _1} = \frac{{ - \sin \mu \sin \left( {\theta + \mu } \right)}}{{\cos \theta }},~~{\lambda _2} = \frac{{ - \sin \mu \cos \left( {\theta + \mu } \right)}}{{\sin \theta }}\\ {\rm{II}}:\;\;\;{\lambda _1} = \frac{{\sin \mu \sin \left( {\theta - \mu } \right)}}{{\cos \theta }},~~{\lambda _2} = \frac{{ - \sin \mu \cos \left( {\theta - \mu } \right)}}{{\sin \theta }}\\ {\rm{III}}:\;\;{\lambda _1} \to \infty ,\;\;\;{\lambda _2} \to \infty \end{array}$

(9)

将这3组解代回式(7)和式(8)，最终得到3组方程：

$\left\{ \begin{array}{l} \frac{{{\rm{d}}y}}{{{\rm{d}}x}}\left| {_{\rm{I}}} \right. = \tan \left( {\theta - \mu } \right)\\ - \frac{1}{{\rho {u^2}}}\frac{{{\rm{d}}p}}{{{\rm{d}}x}}\left| {_{\rm{I}}} \right. = - \tan \mu \frac{{{\rm{d}}\theta }}{{{\rm{d}}x}}\left| {_{\rm{I}}} \right. - \frac{{\cos \theta }}{{\rho {u^2}\cos \mu \cos \left( {\theta - \mu } \right)}}{\left( {\frac{{\partial p}}{{\partial e}}} \right)_\rho }\frac{{T{\rm{d}}s}}{{{\rm{d}}x}}\left| {_{{\rm{III}}}} \right. + \frac{{\delta \sin \theta {{\sin }^2}\mu }}{{y\cos \left( {\theta - \mu } \right)\cos \mu }} \end{array} \right.$

(10)

$\left\{ \begin{array}{l} \frac{{{\rm{d}}y}}{{{\rm{d}}x}}\left| {_{{\rm{II}}}} \right. = \tan \left( {\theta + \mu } \right)\\ - \frac{1}{{\rho {u^2}}}\frac{{{\rm{d}}p}}{{{\rm{d}}x}}\left| {_{{\rm{II}}}} \right. = \tan \mu \frac{{{\rm{d}}\theta }}{{{\rm{d}}x}}\left| {_{{\rm{II}}}} \right. - \frac{{\cos \theta }}{{\rho {u^2}\cos \mu \cos \left( {\theta - \mu } \right)}}{\left( {\frac{{\partial p}}{{\partial e}}} \right)_\rho }\frac{{T{\rm{d}}s}}{{{\rm{d}}x}}\left| {_{{\rm{III}}}} \right. + \frac{{\delta \sin \theta {{\sin }^2}\mu }}{{y\cos \left( {\theta - \mu } \right)\cos \mu }} \end{array} \right.$

(11)

$\left\{ \begin{array}{l} \frac{{{\rm{d}}y}}{{{\rm{d}}x}}\left| {_{{\rm{III}}}} \right. = \tan \theta \\ - \frac{1}{{\rho {u^2}}}\frac{{{\rm{d}}p}}{{{\rm{d}}x}}\left| {_{{\rm{III}}}} \right. = \frac{{{\rm{d}}u}}{{u{\rm{d}}x}}\left| {_{{\rm{III}}}} \right. \end{array} \right.$

(12)

式中：下标“Ⅰ”“Ⅱ”“Ⅲ”分别对应3个求导方向。式(10)代表左行马赫线和线上的相容方程，式(11)代表右行马赫线和线上的相容方程，式(12)代表流线和沿流线的欧拉方程。可以看出，相容方程中的 $\frac{\delta \sin \theta {{\sin }^{2}}\mu }{y\cos \left( \theta \pm \mu \right)\cos \mu }$ 代表轴向变形对流动的影响，而 ${{\left( \frac{\partial p}{\partial e} \right)}_{\rho }}\frac{T\text{d}s}{\text{d}x}\left| _{Ⅲ} \right.$ 为沿流线的熵变影响项，反映熵变对压力的修正。如前所述，由于熵变项分离自压力的随体导数，因此ds/dx所代表的只可能是沿质点轨迹(即流线)的导数。

值得注意的是，在轴线上有y→0, θ→0，因此 $\frac{\delta \sin \theta {{\sin }^{2}}\mu }{y\cos \left( \theta +\mu \right)\cos \ \mu }$ 是“ $\frac{0}{0}$ ”型极限。当y→0时，轴线附近流动可以看作一维变截面管道等熵流，流速u与“管道半径”y有如下关系：

$\frac{{{\rm{d}}u}}{u} = - \frac{1}{{1 - M{a^2}}}\frac{{{\rm{d}}\left( {{\rm{ \mathsf{ π} }}{y^2}} \right)}}{{{\rm{ \mathsf{ π} }}{y^2}}} = \frac{{2{{\tan }^2}\mu {\rm{d}}y}}{y}$

(13)

式中：流体马赫数 $Ma=\frac{u}{c}$ ，将式(11)以及流线方程 $\frac{\text{d}y}{\text{d}x}=\tan \ \theta$ 代入 $\frac{\delta \sin \theta {{\sin }^{2}}\mu }{y\cos \left( \theta \pm \mu \right)\cos \ \mu }$ 中，可得：

$\mathop {\lim }\limits_{\begin{array}{*{20}{c}} {y \to 0}\\ {\theta \to 0} \end{array}} \frac{{\delta \sin \theta {{\sin }^2}\mu }}{{y\cos \left( {\theta \pm \mu } \right)\cos \mu }} = \mathop {\lim }\limits_{\begin{array}{*{20}{c}} {y \to 0}\\ {\theta \to 0} \end{array}} \left( {\frac{{\delta \cos \theta \cos \mu }}{{\cos \left( {\theta \pm \mu } \right)}}\frac{{{{\tan }^2}\mu {\rm{d}}y}}{{y{\rm{d}}x}}} \right) = \frac{{\delta {\rm{d}}u}}{{2u~{\rm{d}}x}}\mathop {\lim }\limits_{\begin{array}{*{20}{c}} {y \to 0}\\ {\theta \to 0} \end{array}} \frac{{\cos \theta \cos \mu }}{{\cos \left( {\theta \pm \mu } \right)}} = \frac{{{\rm{d}}u}}{{2u~{\rm{d}}x}}$

(14)

即该极限趋近于一个有限值，以往做法是取轴线附近一点的值作为替代^[21-22]，实际上取轴线上的du/2u dx即可。

可以看出，式(10)~(12)所示的常微分方程组构成了二维定常可压缩无黏非等熵流的完备三族特征线方程组，在形式上与等熵流特征线方程组^[21]完全一致。其中非等熵的影响体现在相容方程包含的沿流线的熵变项中，相比于以往推导的非等熵特征线方程组^[24]中所指的沿特征线的熵变项，显然沿流线的熵变项才是非等熵流的物理实质。通过合理描述该熵变项，可求解带化学反应或热交换的非等熵流动，也可求解含冲击波的非均熵流(沿流线等熵而流线间存在熵差)，但后者仍需第三族特征线(流线)确定温度、比焓等热力学量。只有对于沿流线没有熵变且流线之间没有熵差的均熵流，三族特征线才能退化为两族特征线方程求解。

2. 柱状装药的水下爆炸模型

考虑柱状装药水下爆炸过程，设装药无限长，爆轰稳定，爆速恒定不变，CJ(Chapman-Jouguet)面垂直于轴线的平面，忽略热传递和界面失稳。若将参考系取在CJ面上，则可建立如图 1所示的定常模型^[25]。爆轰产物的流动是二维定常可压缩均熵无黏柱状流，其中CJ面的马赫数为1^[26]。水中流动是二维定常可压缩非均熵无黏柱状流，水中冲击波是来流速度为爆速的驻定曲冲击波，冲击波强度随着径向距离的增加而衰减，波前后产生熵跳跃，之后沿流线等熵，流线之间存在熵差，为非均熵流。在爆轰产物与水流的水-气界面上，满足压力连续和法向速度连续的边界条件，后者即为水、气两侧偏转角θ相等。在紧邻爆轰波的水-气界面处是曲冲击波的起点，爆轰产物一侧可以用普朗特-迈耶尔(Prandtl-Meyer)绕流描述^[27]：

${\rm{d}}\theta = \sqrt {M{a^2} - 1} \frac{{{\rm{d}}u}}{u}$

(15)

$\rho u~{\rm{d}}u + {\rm{d}}p = 0$

(16)

图 1 柱状装药水下爆炸模型

Figure 1. Underwater explosion model of cylindrical charge

下载: 全尺寸图片幻灯片

式中：Ma=u/c为爆轰产物的马赫数。由于爆轰产物为均熵流，因此仅在流线上满足的式(12)转化为全流场适用的式(16)。再将水中斜冲击波关系与爆轰产物普朗特-迈耶尔绕流方程联立，即可求出流场的初始值。

具体计算时，对炸药的爆轰产物选用JWL(Jones-Wilkins-Lee)状态方程描述^[28]：

$p = A\left( {1 - \frac{\omega }{{{R_1}V}}} \right){{\rm{e}}^{ - {R_1}V}} + B\left( {1 - \frac{\omega }{{{R_2}V}}} \right){{\rm{e}}^{ - {R_2}V}} + \frac{{\omega E}}{V}$

(17)

式中：v₀为炸药初始比容，V=v/v₀为相对体积；E为体积内能；A、B、R₁、R₂、ω为常数，且(∂E/∂p)_ρ=V/ω。对水选用如下形式的Mie-Grüneison状态方程描述：

$p = \left\{ \begin{array}{l} {A_1}\mu + {A_2}{\mu ^2} + {A_3}{\mu ^3} + \left( {{B_0} + {B_1}\mu } \right){\rho _0}e\;\;\;\;\;\;\;\;\mu \ge 0\\ {T_1}\mu + {T_2}{\mu ^2} + {B_0}{\rho _0}e\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mu < 0 \end{array} \right.$

(18)

式中：ρ₀为水的初始密度，μ为压缩率，μ=ρ/ρ₀-1；e为比内能，A₁、A₂、A₃、T₁、T₂、B₀、B₁均为系数，且有：

${\left( {\frac{{\partial e}}{{\partial p}}} \right)_\rho } = \left\{ \begin{array}{l} 1/\left[ {\left( {{B_0} + {B_1}\mu } \right){\rho _0}} \right]\;\;\;\;\mu \ge 0\\ 1/\left( {{B_0}{\rho _0}} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mu < 0 \end{array} \right.$

(19)

状态方程的具体参数见表 1和表 2，其中炸药PETN、TNT以及水的参数源于文献[29]，SEP炸药(其中PETN和石蜡的质量分数分别为65%和35%)的参数源于文献[4]，D为爆速，p_CJ为CJ压力。在等熵卸载过程中，爆轰产物或水需要使用等熵条件：

$T{\rm{d}}s\left| {_{{\rm{III}}}} \right. = \left( {{\rm{d}}e + p{\rm{d}}v} \right)\left| {_{{\rm{III}}}} \right. = 0$

(20)

表 1 几种炸药的爆轰产物JWL状态方程参数

Table 1. Parameters of JWL equation of state for several explosives

炸药	CJ参数			JWL状态方程参数
炸药	ρ₀/(g·cm^-3)	D/(km·s^-1)	p_CJ/GPa	A/GPa	B/GPa	R₁	R₂	ω
PETN 0.88	0.880	5.170	6.2	348.62	11.29	7.0	2.0	0.24
PETN 1.26	1.260	6.540	14.0	573.10	20.16	6.0	1.8	0.28
PETN 1.50	1.500	7.450	22.0	625.30	23.29	5.25	1.6	0.28
PETN 1.77	1.770	8.300	33.5	617.05	16.93	4.4	1.2	0.25
TNT	1.630	6.930	21.0	373.77	3.75	4.15	0.9	0.35
SEP	1.310	6.970	15.9	372.00	3.48	4.59	1.06	0.29

下载: 导出CSV

| 显示表格

表 2 水的Mie-Grüneison状态方程

Table 2. Mie-Grüneison equation of state of water

ρ₀/(g·cm^-3)	e₀/(J·kg^-3)	p_CJ/GPa	A₁/GPa	A₂/GPa	A₃/GPa	B₀	B₁	T₁/GPa	T₂/GPa
1.000	361.9	6.2	2.2	9.54	14.57	0.28	0.28	2.2	0

下载: 导出CSV

| 显示表格

3. 水下爆炸计算实例

3.1 计算格式

非等熵流的三族特征线方程组实际上就是非齐次常微分方程组。以4节点计算格式(见图 2(a))为例，计算格式为：

$\frac{{y_j^{i + 1} - y_{j \pm 1}^i}}{{x_j^{i + 1} - x_{j \pm 1}^i}} = {\left[ {\tan \left( {\theta \pm \mu } \right)} \right]_{{\rm{I}},{\rm{II}}}}$

(21)

$\begin{array}{l} {\left[ { - \frac{1}{{\rho {u^2}}}} \right]_{{\rm{I}},{\rm{II}}}}\frac{{p_j^{i + 1} - p_{j \pm 1}^i}}{{x_j^{i + 1} - x_{j \pm 1}^i}} = {\left[ { \pm \tan \mu } \right]_{{\rm{I}},{\rm{II}}}}\frac{{\theta _j^{i + 1} - \theta _{j \pm 1}^i}}{{x_j^{i + 1} - x_{j \pm 1}^i}} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left[ {\frac{{T\cos \theta }}{{\rho {u^2}\cos \mu \cos \left( {\theta \pm \mu } \right)}}{{\left( {\frac{{\partial p}}{{\partial e}}} \right)}_\rho }} \right]_{{\rm{III}}}}\frac{{s_j^{i + 1} - s_j^i}}{{x_j^{i + 1} - x_j^i}} + {\left[ {\frac{{\delta \sin \theta {{\sin }^2}\mu }}{{y\cos \left( {\theta \pm \mu } \right)\cos \mu }}} \right]_{{\rm{I}},{\rm{II}}}} \end{array}$

(22)

图 2 二维特征线差分法的计算格式

Figure 2. Computational scheme of two-dimensional FDM of characteristics

下载: 全尺寸图片幻灯片

式中：“[]_{Ⅰ, Ⅱ}”代表沿特征线的加权取值；i、j为索引序号，代表逻辑空间的网格。为了避免隐式计算，常用显格式的欧拉预估-校正法(经典特征线理论称参数平均法^[30])，即第一步预估的“[]_{Ⅰ, Ⅱ}”取已知点的参数计算解点，第二步校正的“[]_{Ⅰ, Ⅱ}”取已知点和解点参数的平均值，预估-校正法在形式上具有二阶精度(Δx², Δy²)。更普遍的5节点计算格式(见图 5(b))用4节点的计算格式通过内插获得。

图 5 柱状TNT水下爆炸在R=2R₀处的压力时程曲线

Figure 5. Pressure history curve at R=2R₀ in underwater explosion of cylindrical TNT

下载: 全尺寸图片幻灯片

3.2 计算结果

在特征线差分法中，曲冲击波的形状通过一系列斜冲击波逼近，逼近程度与差分步长相关。冲击波形状是压力、速度等发生突跃变化的不连续面，用有限元软件如AUTODYN也能大致观测到。图 3展示了4种炸药的特征线差分法计算结果(Cal., PETN)和AUTODYN模拟结果(计算区域1 500 mm×600 mm，网格尺寸1 mm×1 mm)，其中R、x为径向和轴向距离，R₀为装药半径。图 4展示了SEP炸药的特征线差分法计算结果和实验结果^[4]。可以看出：特征线差分法比较准确地捕捉到了近场冲击波形状。

图 3 4种柱状PETN炸药水下爆炸冲击波形状

Figure 3. Underwater shock wave shape of four kinds of cylindrical explosives PETN

下载: 全尺寸图片幻灯片

图 4 柱状SEP炸药水下爆炸冲击波形状

Figure 4. Underwater shock wave shape of cylindrical explosive SEP

下载: 全尺寸图片幻灯片

根据斜冲击波理论^[21]，斜冲击波上的物理参数如压力、速度、偏转角、激波角等，只有一个参数是独立的。如果已知曲冲击波形状，那么从数学上可以确定曲冲击波面上每一处的斜率(即激波角)，从而确定整个曲冲击波面上的压力、速度、偏转角等。如果特征线差分法计算的近场冲击波形状足够精确，可以推断：特征线差分法计算的近场冲击波压力也足够精确。图 5展示了特征线差分法和AUTODYN所代表的有限元法计算的某空间位置(R=2R₀)的压力时程曲线。可以看出：有限元法得到的压力时程曲线具有明显波动，而特征线差分法获得的压力时程曲线可以瞬间跳跃、连续衰减而无任何波动。

4. 结论

根据二维定常可压缩超声速非等熵流的控制方程，推导出其特征线方程，并通过补充流线方程作为第三族特征线方程，使特征线差分方法可以求解非等熵流问题，相比于以往添加沿特征线的熵变项，更能体现非等熵流的物理实质。

采用三族特征线方法，将非等熵流问题转化为求解非齐次常微分方程组问题，提出了五点差分格式及其欧拉预估-校正解法，在理论上可以保证二阶计算精度。

对柱状装药的水下爆炸建立了定常模型，将非等熵特征线差分法应用于求解水下爆炸近场的非均熵流动。对几种炸药的水下爆炸近场冲击波进行计算，结果表明，新的特征线差分方法可以准确地捕捉冲击波形状并获得波后冲击波压力历程，避免了常规方法计算冲击波压力的人工黏性误差和数值波动，说明所提出的三族特征线差分法可以用于处理水下爆炸这种非均熵问题。

图 1 柱状装药水下爆炸模型

Figure 1. Underwater explosion model of cylindrical charge

下载: 全尺寸图片幻灯片

图 2 二维特征线差分法的计算格式

Figure 2. Computational scheme of two-dimensional FDM of characteristics

下载: 全尺寸图片幻灯片

图 5 柱状TNT水下爆炸在R=2R₀处的压力时程曲线

Figure 5. Pressure history curve at R=2R₀ in underwater explosion of cylindrical TNT

下载: 全尺寸图片幻灯片

图 3 4种柱状PETN炸药水下爆炸冲击波形状

Figure 3. Underwater shock wave shape of four kinds of cylindrical explosives PETN

下载: 全尺寸图片幻灯片

图 4 柱状SEP炸药水下爆炸冲击波形状

Figure 4. Underwater shock wave shape of cylindrical explosive SEP

下载: 全尺寸图片幻灯片

表 1 几种炸药的爆轰产物JWL状态方程参数

Table 1. Parameters of JWL equation of state for several explosives

炸药	CJ参数			JWL状态方程参数
炸药	ρ₀/(g·cm^-3)	D/(km·s^-1)	p_CJ/GPa	A/GPa	B/GPa	R₁	R₂	ω
PETN 0.88	0.880	5.170	6.2	348.62	11.29	7.0	2.0	0.24
PETN 1.26	1.260	6.540	14.0	573.10	20.16	6.0	1.8	0.28
PETN 1.50	1.500	7.450	22.0	625.30	23.29	5.25	1.6	0.28
PETN 1.77	1.770	8.300	33.5	617.05	16.93	4.4	1.2	0.25
TNT	1.630	6.930	21.0	373.77	3.75	4.15	0.9	0.35
SEP	1.310	6.970	15.9	372.00	3.48	4.59	1.06	0.29

下载: 导出CSV

表 2 水的Mie-Grüneison状态方程

Table 2. Mie-Grüneison equation of state of water

ρ₀/(g·cm^-3)	e₀/(J·kg^-3)	p_CJ/GPa	A₁/GPa	A₂/GPa	A₃/GPa	B₀	B₁	T₁/GPa	T₂/GPa
1.000	361.9	6.2	2.2	9.54	14.57	0.28	0.28	2.2	0

下载: 导出CSV

参考文献(30)

[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. https://www.researchgate.net/publication/236565135_Spherical_explosions_and_the_equation_of_state_of_water
[3]	池家春, 马冰.TNT/RDX(40/60)炸药球水中爆炸波研究[J].高压物理学报, 1999, 13(3):199-204. doi: 10.11858/gywlxb.1999.03.008 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. doi: 10.11858/gywlxb.1999.03.008
[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. http://www.researchgate.net/publication/241230429_Investigations_of_fundamental_properties_of_underwater_shock_waves_by_high-speed_photography
[5]	赵继波, 谭多望, 李金河, 等.TNT药柱水中爆炸近场压力轴向衰减规律[J].爆炸与冲击, 2008, 28(6):539-543. doi: 10.11883/1001-1455(2008)06-0539-05 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. doi: 10.11883/1001-1455(2008)06-0539-05
[6]	张振华, 朱锡, 白雪飞.水下爆炸冲击波的数值模拟研究[J].爆炸与冲击, 2004, 24(2):182-188. http://www.bzycj.cn/CN/abstract/abstract9940.shtml 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. http://www.bzycj.cn/CN/abstract/abstract9940.shtml
[7]	方斌, 朱锡, 张振华, 等.水下爆炸冲击波数值模拟中的参数影响[J].哈尔滨工程大学学报, 2005, 26(4):419-424. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=hebgcdxxb200504001 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. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=hebgcdxxb200504001
[8]	胡毅亭, 贾宪振, 饶国宁, 等.水下爆炸冲击波和气泡脉动的数值模拟研究[J].舰船科学技术, 2009, 31(2):134-140. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=jckxjs200902032 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. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=jckxjs200902032
[9]	刘科种, 徐更光, 辛春亮, 等.AUTODYN水下爆炸数值模拟研究[J].爆破, 2009, 26(3):18-21. http://www.cqvip.com/qk/96961X/200903/31592781.html LIU Kezhong, XU Gengguang, XIN Chunliang, et al. Research on numerical simulation in underwater explosion by AUTODYN[J]. Blasting, 2009, 26(3):18-21. http://www.cqvip.com/qk/96961X/200903/31592781.html
[10]	辛春亮, 秦健, 刘科种, 等.基于LS-DYNA软件的水下爆炸数值模拟研究[J].弹箭与制导学报, 2008, 28(3):156-158. http://www.cqvip.com/QK/97591A/200803/27650534.html 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. http://www.cqvip.com/QK/97591A/200803/27650534.html
[11]	汪俊, 刘建湖, 李玉节.加筋圆柱壳水下爆炸动响应数值模拟[J].船舶力学, 2006, 10(2):126-137. http://cdmd.cnki.com.cn/Article/CDMD-10487-2009033857.htm 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. http://cdmd.cnki.com.cn/Article/CDMD-10487-2009033857.htm
[12]	JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126(1):202-228. doi: 10.1006/jcph.1996.0130
[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. doi: 10.1016/j.jcp.2010.08.016
[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. doi: 10.1007/s00466-002-0371-6
[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. http://www.sciencedirect.com/science/article/pii/S004579301200401X
[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. doi: 10.1006/jcph.1994.1155
[17]	师华强, 汪玉, 宗智, 等.近水面水下爆炸二维Level-set数值模拟[J].计算力学学报, 2010, 27(3):409-414. doi: 10.7511/jslx20103005 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. doi: 10.7511/jslx20103005
[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. doi: 10.2514/3.4038
[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. http://www.researchgate.net/publication/23926144_Comparison_of_characteristics_and_shock_capturing_methods_with_application_to_the_space_shuttle_vehicle
[22]	ZUCROW M J, HOFFMAN J D. 气体动力学(上册)[M]. 王汝涌, 吴宗真, 吴宗善, 等译. 北京: 国防工业出版社, 1984: 431-432.
[23]	李晓杰, 张程娇, 闫鸿浩, 等.水下爆炸近场非均熵流的特征线差分解法[J].爆炸与冲击, 2012, 32(6):604-608. doi: 10.11883/1001-1455(2012)06-0604-05 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. doi: 10.11883/1001-1455(2012)06-0604-05
[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. https://www.researchgate.net/publication/236480548_Some_practical_methods_of_using_characteristics_in_the_calculation_of_non-steady_compressible_flow
[25]	STEBNOVSKⅡ 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. doi: 10.1007/BF00783052
[26]	李晓杰, 王金相, 闫鸿浩, 等.基于Laval喷管的圆管爆轰驱动近似解[J].工程爆破, 2003, 9(4):7-9. http://cpfd.cnki.com.cn/Article/CPFDTOTAL-GCWL199800001021.htm 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. http://cpfd.cnki.com.cn/Article/CPFDTOTAL-GCWL199800001021.htm
[27]	李晓杰, 赵春风.基于通用炸药状态方程分析飞板运动规律的特征线法[J].爆炸与冲击, 2012, 32(3):237-242. doi: 10.11883/1001-1455(2012)03-0237-06 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. doi: 10.11883/1001-1455(2012)03-0237-06
[28]	周清, 张云海, 毛允波.典型民用建筑内爆炸波传播机理[J].山东科技大学学报(自然科学版), 2009, 28(6):36-44. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=sdkjdxxb200906007 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. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=sdkjdxxb200906007
[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.

施引文献

期刊类型引用(2)

1.	李晓杰，杨晨琛，闫鸿浩，王小红，王宇新，张程娇. 柱状含铝炸药水下爆炸近场的特征线法研究. 爆炸与冲击. 2019(02): 22-29 . 本站查看
2.	杨晨琛，李晓杰，闫鸿浩，王小红，王宇新. 爆轰产物状态方程的水下爆炸反演理论研究. 爆炸与冲击. 2019(09): 26-36 . 本站查看

其他类型引用(1)

资源附件(0)

访问统计

图(5) / 表(2)

计量

文章访问数: 5739
HTML全文浏览量: 2168
PDF下载量: 245
被引次数: 3

1. 二维定常非等熵流的特征线方程
2. 柱状装药的水下爆炸模型
3. 水下爆炸计算实例
3.1 计算格式
3.2 计算结果
4. 结论

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

doi: 10.11883/bzycj-2016-0314

作者简介: 李晓杰(1963-), 男, 教授, 博士生导师, dalian03@qq.com

计量

出版历程