2024年5月26日发(作者：仇英武)
第３３卷第６期
２００９年ｌ２月
南京理工大学学报（自然科学版）
Ｊｏｕｒｎａｌ ｏｆ Ｎａｎｊｉｎｇ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｓｃｉｅｎｃｅ ａｎｄ Ｔｅｃｈｎｏｌｏｇｙ（Ｎａｔｕｒａｌ Ｓｃｉｅｎｃｅ）
Ｖ０１．３３ Ｎｏ．６
Ｄｅｃ．２００９
冲击波超压峰值的数值计算
王 杨，郭则庆，姜孝海
（南京理工大学瞬态物理国家重点实验室，江苏南京２１００９４）
摘要：为提高冲击波超压峰值的数值计算精度，该文根据二阶Ｒｏｅ格式，模拟了３种不同粗细
网格条件下的点爆炸流场。基于高精度的数值计算格式，利用冲击波的发展规律求得波阵面速
度，再由兰金一雨贡纽方程计算得到超压峰值。并与数值模拟的直接超压峰值进行了对照讨论。
结果表明随着网格加密，数值直接超压峰值不断逼近采用冲击波传播速度获得的超压峰值．而
后者几乎不变。
关键词：冲击波；Ｒｏｅ格式；冲击波阵面；兰金一雨贡纽方程；超压峰值
中图分类号：ＴＪ ０１２．２ 文章编号：１００５—９８３０（２００９）０６—０７７０—０４
Ｎｕｍｅｒｉｃａｌ Ｉｎｖｅｓｔｉｇａｔｉｏｎｓ ｏｎ Ｐｅａｋ Ｏｖｅｒｐｒｅｓｓｕｒｅ ｏｆ Ｂｌａｓｔ Ｗａｖｅ
ＷＡＮＧ Ｙａｎｇ，ＧＵＯ Ｚｅ—ｑｉｎｇ，ＪＩＡＮＧ Ｘｉａｏ—ｈａｉ
（Ｎａｔｉｏｎａｌ Ｋｅｙ Ｌａｂｏｒａｔｏｒｙ ｏｆ Ｔｒａｎｓｉｅｎｔ Ｐｈｙｓｉｃｓ，ＮＵＳＴ，Ｎａｎｊｉｎｇ ２１００９４，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ｉｎ ｏｒｄｅｒ ｔｏ ｉｍｐｒｏｖｅ ｔｈｅ ｃｏｍｐｕｔａｔｉｏｎ ｐｒｅｃｉｓｉｏｎ ｏｆ ｐｅａｋ ｏｖｅｒｐｒｅｓｓｕｒｅ ｏｆ ｔｈｅ ｂｌａｓｔ ｗａｖｅ，
ｔｈｒｅｅ ｐｏｉｎｔ ｅｘｐｌｏｓｉｏｎ ｃａｓｅｓ ａｒｅ ｓｉｍｕｌａｔｅｄ ｎｕｍｅｒｉｃａｌｌｙ ｗｉｔｈ ｄｉｆｆｅｒｅｎｔ ｆｉｎｅ ｇｒｉｄｓ ｂｙ ｕｓｉｎｇ ｔｈｅ ｓｅｃｏｎｄ－－ｏｒ－・
ｄｅｒ Ｒｏｅ ｓｃｈｅｍｅ．Ｂａｓｅｄ ｏｎ ｔｈｅ ｎｕｍｅｒｉｃａｌ ｃｏｍｐｕｔａｔｉｏｎａｌ ｓｃｈｅｍｅ ｏｆ ｔｈｅ ｈｉｇｈ—ｐｒｅｃｉｓｉｏｎ，ｔｈｅ ｂｌａｓｔ
ｐｒｏｐａｇａｔｉ’ｏｎ ｓｐｅｅｄ ｉｓ ｃａｌｃｕｌａｔｅｄ ａｃｃｏｒｄｉｎｇ ｔｏ ｔｈｅ ｄｅｖｅｌｏｐｍｅｎｔ ｌａｗ ｏｆ ａ ｂｌａｓｔ ｗａｖｅ．Ｔｈｅ ｉｎｄｉｒｅｃｔ ｐｅａｋ
ｏｖｅｒｐｒｅｓｓｕｒｅ ｉｓ ｃａｌｃｕｌａｔｅｄ ｂｙ ｕｓｉｎｇ Ｒａｎｋｉｎｅ－－Ｈｕｇｏｎｉｏｔ ｅｑｕａｔｉｏｎ ａｎｄ ｃｏｍｐａｒｅｄ ｗｉｔｈ ｔｈｅ ｄｉｒｅｃｔ ｐｅａｋ Ｏ—－
ｖｅｒｐｒｅｓｓｕｒｅ ｃａｌｃｕｌａｔｅｄ ｂｙ ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ．Ｔｈｅ ｃｏｍｐａｒａｔｉｖｅ ｒｅｓｕｌｔ ｉｎｄｉｃａｔｅｓ ｔｈａｔ：ｗｉｔｈ ｔｈｅ ｉｎ。
ｃｒｅａｓｅ ｏｆ ｔｈｅ ｍｅｓｈ ｒｅｆｉｎｅｍｅｎｔ，ｔｈｅ ｄｉｒｅｃｔ ｎｕｍｅｒｉｃａｌ ｐｅａｋ ｏｖｅｒｐｒｅｓｓｕｒｅ ｇｒａｄｕａｌｌｙ ａｐｐｒｏａｃｈｅｓ ｔｈｅ ｉｎｄｉ－
ｒｅｃｔ ｐｅａｋ ｏｖｅｒｐｒｅｓｓｕｒｅ ｏｂｔａｉｎｅｄ ｂｙ ａ ｂｌａｓｔ ｐｒｏｐａｇａｔｉｏｎ ｓｐｅｅｄ，ｗｈｉｌｅ ｔｈｅ ｌａｔｔｅｒ ｉｓ ａｌｍｏｓｔ ｃｏｎｓｔａｎｔ．
Ｋｅｙ ｗｏｒｄｓ：ｂｌａｓｔ ｗａｖｅｓ；Ｒｏｅ ｓｃｈｅｍｅ；ｓｈｏｃｋ ｆｒｏｎｔｓ；Ｒａｎｋｉｎｅ－Ｈｕｇｏｎｉｏｔ ｅｑｕａｔｉｏｎ；ｐｅａｋ ｏｖｅｒｐｒｅｓ‘
Ｓｌ１ｒＰ
火炮发射时，高温、高压火药燃气在弹丸飞出
冲击波超压峰值的测量精度与计算精度就成了本
后，自膛Ｅｌ急剧膨胀，形成膛口冲击波，对人员、装
备及周围设施形成了严重的危害。为此，各国军方
都制定了相应的安全和防护标准¨＇２ Ｊ，作为武器定
领域一直在探讨的理论与技术课题 。由于冲
击波阵面是一个间断面，因此超压峰值的测量与计
算精度受到测量系统响应与计算格式及网格的影
型与装备使用的依据之一。而冲击波的超压峰值
与冲量是该防护标准的２个基本参量。于是提高
收稿日期：２００９—０８—１０ 修回日期：２００９—１０—２７
响。近年来，一直在寻求提高精度的方法。
在试验测量方面，因冲击波超压测试系统带
基金项目：国防科技重点实验室基金
作者简介：王杨（１９７２一），女，博士生，主要研究方向：偏微分方程数值解，兵器发射理论与技术，Ｅ－ｍａｉｌ：ｗａｎｇ），一
ａｎｇｘｕｐｕ＠１６３．ｃｏｎ。
总第１６９期 王杨郭则庆姜孝海 冲击波超压峰值的数值计算 ７７１
宽有限，使得由信号响应引起的测量误差为５％
～
１０％，而这种误差一般与冲击波强度直接相关，
因此，系统实时修正及现场标定就成了保证冲击
波超压峰值测量精度的基本途径？在冲击波测量
系统标定方法的研究中，人们提出了不同的标定
方法：文献［６—８］针对不同类型的传感器提出
了其实验室标定方法和现场标定方法：美国陆军
试验规程 巾规定了利用冲击波速度标定超压
峰值的方法
对于满足物理方程的简单模型（例如不考虑
外界反应的强爆炸冲击波模型），可以得到解析
解，其超压峰值的计算是精确的 但是，描述复杂
流场冲击波（如考虑反压的弱冲击波和反射、绕
射冲击波等）运动规律的控制方程是非线性偏微
分方程组，一般多采用数值计算进行离散求解：
只有在极少情况下，事先预知激波位置，可以用
“激波装配法”得到准确的超压峰值。多数情况
下均采用“激波捕捉法”，其数值解的精度与所选
取的格式及网格粗细有关：而计算格式中加入的
人 １ 粘性，在数值求解中抹平了间断 ：选取间
断面附近不同的格点，将会得到不同的超压峰值
（文中将这种由计算格点直接获取的超压峰值称
为直接超压峰值，而由数值结果通过某种方法间
接获得的超压峰值，称为问接超压峰值） 在一
定意义上，随着网格间距的减小，直接超压峰值的
精度会有所提高。然而，随着冲击波的传播（尤
其是冲击波的远场传播），其计算区域迅速扩大，
计算量剧增，对计算机的计算速度、存储量等提出
了更高的要求；虽然存在局部加密方法（如自适
应网格等），但计算量与计算时问的矛盾仍比较
突出：如何采用较少的网格，计算 足够精度的
超压峰值就成为迫切需要解决的问题：
本文基于高精度的数值计算格式，利用冲击
波的发展规律求得波阵面速度，再由兰金一雨贡纽
方程计算得到超压峰值。根据不同网格疏密条件
下的计算结果，讨论了直接超压峰值和间接超压
峰值的变化特点＝
基本方程和计算格式
在不考虑粘性和化学反应的条件下，理想气
体满足
杀 Ｑｄ
¨ｔ）
Ｖ＋
ｆｔ）
ｄＳ＝０ （１）
式中：
Ｑ：［Ｐ ｐＵ Ｅｒ
Ｆ＝Ｑ ＋Ｇ
Ｇ＝［０ Ｐｎ ＰＵ］
ｐ、Ｕ、Ｅ分别表示流体密度、速度和单位体积的内
能，Ｓ（ｔ）表示ｔ时刻控制体Ｖ（ｔ）的表面积，ｎ表示
其外法线方向。为使方程闭合，引入状态方程
Ｐ＝ＪＤ（ —１）ｅ
式中：比热比ｙ＝１．４，ｅ为比内能。
将方程（１）离散，得
Ａ ｔ
＋
主 Ｓ：
＿一
‘‘
ｉ＝
Ｉ
０
即
＝
（ 古 Ｎ ＦｉＳ Ａ ） （２）
式中：上标ｎ表示前一时刻，ｎ＋１表示所求的时
刻。Ⅳ表示所求控制体所包含的外表面总数，Ｓ
表示第ｉ个面的面积，Ｆ 表示Ｆ在该面外法线方
向的分量，即界面流量；将控制体界面视为一维
黎曼问题，界面流量根据Ｒｏｅ方法求解
Ｆ ＝Ｆ （ＱＲ，ＱＩ）Ａ＝ＯＦ／ＯＱ
Ｆ ＝［Ｆ（Ｑ ）＋Ｆ（Ｑ ）一∑ ＩＡ，Ｉ ｅ，］／２
式中：Ｑ 和Ｑ 分别表示界面左右两侧格点值，Ａ
为Ｅｕｌｅｒ方程的Ｊａｃｏｂｉａｎ矩阵， ，Ａ ，ｅｊ分别表示
波强度矩阵 、特征向量矩阵 和右特征矩阵日
的分量。
为提高精度，引入二阶修正项
：
∑ （１一 Ａｔ ）ｗ ／２
式中：
Ｑ 一Ｑ ＝∑Ｗｉ ∑ｏｌｉ
因此，界面流量为
Ｆ ＝［Ｆ（Ｑ ）＋Ｆ（Ｑ ）一∑ ＩＡ ｌｅｊ］／２＋
∑ （１一 Ａｔ ）Ｗｉ／２
在跨声速或超声速计算中，当流场产生间断
时（如激波），其周围会产生所谓的伪震荡（Ｓｐｕｒｉ—
ＯｔＩＳ ｏｓｃｉｌｌａｔｉｏｎ） ，出现非物理解。因而需要对
其进行修正，即
Ｗｉ＝ Ｉｒ
式中：
＝（ｂ（０ ）
Ａ 如
：
Ｉ：
Ｏｔｉ Ｉ ＋ｌ Ａ ≤０
７７２ 南京理工大学学报（自然科学版） 第３３卷第６期
关于 （０）的取值，本文采用Ｍｏｎｏｔｏｎｉｚｅｄ Ｃｅｎｔｅｒｅｄ
方法，即
（０）＝ｍａｘ（０，ｒａｉｎ（（１＋０）／２，２，２０））
如前文所述，称为间接超压峰值。为了计算超压
峰值，只需在一般流场计算完毕后，再调用计算超
压峰值的子程序进行处理即可。
求得的界面流量具有二阶精度，代入式（２），即可
求得新时刻的值。对于数值解的时间步采用二阶
精度的Ｒｕｎｇｅ．Ｋｕｔｔａ法。
本文考察的点爆炸流场初始条件为：半径为
３算例与讨论
２冲击波超压峰值计算方法
根据前述数值方法获得的流场数值解，利用
波阵面上压强梯度变化最大的特征来确定波阵面
的位置，再由波阵面的位移△ｒ和相应时间间隔
△ ，可求得在此时间间隔内的平均速度。由于
一
０．１ ｍ的静止气球中充满高温、高压的高能气体，
其平均压力为２５６ｐ。，平均密度为８２ｐ。，外部为标
准条件下的大气环境，即Ｐ ＝０．１ ＭＰａ，密度为
Ｐ ＝１．１８ ｋｇ／ｉｎ 。整个计算域为半径１ １ＴＩ的球
体，外边界为外推边界。
分别模拟了网格间距为１／５００ ｍ、１／１ ０００ ＩＴＩ、
般很小（本文Ａｔ＝０．６８２５×１０ Ｓ），因此波阵
１／１ ５００ ｍ的３种情形，对应算例为Ｔ１、 、Ｔ３。
在获得数值解之后，再调用前述的计算超压峰值
的子程序分别计算间接超压峰值Ｐ （见表１）。而
直接超压峰值的取值格点如图１所示，即直接超
压峰值参考Ｐ ，Ｐ 两格点处的超压值。
面速度Ｄ可以表示为
１ ３： ：
ｔ— ●ｄ ｔ ≈——＝—Ａｔ
，一ｔ
二 （３）
、
【ｊ】
式中：ｒ。，ｒ 分别为波阵面在ｔ ，ｔ 时刻的位置。
利用激波关系式
丝：
ｐ１ ｙ＋ｌ
一
ｙ＋ｌ
式中：ｐ 、Ｍａ。＝Ｄ／ｃ 、ｃ。分别为波前气体压力、波阵
面马赫数、声速，Ｐ，为波后气体压力，进一步整理得
＝
（Ｍａ２￣－１）△ｐ＝Ｐ —Ｐ （４）
通过式（４）可以求出冲击波的超压峰值，这种由
波阵面速度和激波关系式间接获取的超压峰值，
图１ 直接超压峰值的选取示意图
表１ 算例Ｔｌ～１、３间接超压峰值与直接超压峰值列表
由表１可以看出，随着网格问距的减小，直接
为了考察直接超压峰值的精度与网格粗细的
关系，计算了算例Ｔ１、Ｔ２、Ｔ３中的直接超压峰值
Ｐ 与算例１、３的间接超压峰值Ｐ 的相对误差Ｅ 、
Ｅ，、Ｅ （见表２），随着网格的加密，相对误差越来
越小，说明直接超压峰值的精度受网格粗细的影
超压峰值Ｐ 、Ｐ 越来越逼近间接超压峰值Ｐ ，而
间接超压峰值并未受网格间距的影响，其值Ｐ 变
化很小，趋于恒定。也就是说在网格较粗的条件
下，通过间接方法得到的超压峰值就能获得较高
的精度（密网格条件下的直接超压峰值）。
响较大。为了得到较精确的直接超压峰值，必须
总第１６９期 王 杨 郭则庆 姜孝海 冲击波超压峰值的数值计算 ７７３
加密网格。但必须注意到，由于加密网格而引起
的计算量、存储量以及计算时间等将会大大增加：
可见，在本文的研究范围内，采用间接方法计算超
压峰值较直接方法（加密网格）更有现实意义和
可行性
表２算例Ｔ１～１、３直接超压峰值与 间接
超压峰值的相对误差列表
图２为算例１、３的在不同时刻的压力和径向
速发典型分布图：从图２（ａ）可见，初始冲击波形
成的同时，膨胀波从接触面开始朝球心内部运动
球心内部点的压力初始保持定值，膨胀波经过后，
压力降低，形成第二道激波，见图２（ｂ）：图２（ｃ）
展示的是典型的压力图，主冲击波后是负相区：
１・４
１ ２
日
皇０．８
运
Ｏ４
０
１５０ ２ｏ０ ２５０ ３Ｏ０ ３５０
ｍｍ
（ａ）ｎ＝１
ｆｂ）”，＝６
３５
２
３０
２５
１
２０ ０
１５
一
１
１０
－
０５
２
６００ ８００
ｒ／ｎｌｒｆｌ
（ｃ）ｎ ＝１６
图２不同时刻的压力与速度分布曲线
（ｎ，表示时间步数）
４结束语
本文采用二阶精度Ｒｏｅ格式，模拟了３种不同
粗细网格的点爆炸流场：利用激波关系式和波阵
面速度间接获得了冲击波的超压峰值 与直接超
压峰值的对比分析表明：采用间接计算方法，可以
在粗网格条件下获得比较准确的初始冲击波超压
峰值 这在一定意义上解决了计算区域大、网格数
量多造成的计算机时长与精度之间的矛盾：
参考文献
『１］ ＭＩＩ 一ＳＴＤ １４７４（ＭＩ）一１９７３，Ｎｏｉｓｅ Ｌｉｍｉｔｓ ｆｏｒ ａｒｍ ｙ ｍａ—
ｔｅｒｉｅｌ［ｓ］．
『２］ ＧＪＢ １１５８—９１，炮口冲击波对人员非听觉器官损伤
的安全限值［ｓ］．
『３ １ Ｋｉｎｎｅｖ Ｇ Ｆ．Ｇｒａｈａｍ Ｋ ｊ．Ｅｘｐｌｏｓｉｖｅ ｓｈｏｃｋｓ ｉｎ ａｉｒ
［Ｍ］．Ｂｅｒｌｉｎ ａｎｄ Ｎｅｗ Ｙｏｒｋ：Ｓｐｒｉｎｇｅｒ—Ｖｅｒｌａｇ，１ ９８５．
［４］ Ｖａｎｄｅｒｓｔｒａｅｔｅｎ Ｂ，ｌ￣ｅｆｅｂｖｒｅ Ｍ，Ｂｅｒｇｈｍａｎｓ Ｊ．Ａ ｓｉｍ‘
ｐｉｅ ｂｌａｓｔ ｗａｖｅ ｍｏｄｅｌ ｆｏｒ ｂｕｒｓｔｉｎｇ ｓｐｈｅｒｅｓ ｂａｓｅｄ Ｏｎ ｎＵ—
ｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ『Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ｈａｚａｒｄｏｕｓ Ｍａｔｅｒｉ—
ａｌｓ，１９９６，４６（２）：１４５—１５７．
［５］李鸿志．炮口制退器科研阶段报告［Ｒ 南京：华东
工程学院，１９７８．１０．
［６］ 张远平，池家春．爆炸冲击波压力传感器灵敏度的
动态标定及测试技术研究［Ａ］．第四届全 爆炸力
学实验技术学术会议论文集［Ｃ ．绵阳：中国力学
学会，２００６．３４９—３５３．
［７］ 崔海涛，刘庆明．冲击波压力传感器测试系统的动
态标定［Ｊ］．流体力学实验与测量，２００４，１８（１）：９２
—
９６．
［８］ 王等旺，张德志，李焰，等．冲击波压扦式传感器测
试系统的动态标定［Ａ］．第十二届全国激波与激
波管学术会议论文集［Ｃ ．洛阳：中国力学学会，
２００６．３６３—３６６．
［９］ ＭＴＰ ４—３—８２２，美陆军试验鉴定部兵器试验规程一
冲击波压力测量（电子的）技术试验［Ｓ ．
ｆ ｌ０ １ Ｍｕｒｄｏｃｋ Ｊ Ｗ．Ｓｈｏｃｋ—ｗａｖｅ ｉｎｔｅｒａｃｔｉｏｎ Ｍｔｈ ｔｗｏ．ｄｉｍｅｎ—
ｓｉｏｎａｌ ｂｏｄｉｅｓ［Ｊ］．ＡＩＡＡ Ｊｏｕｒｎａｌ，１９７５．１３（９），１１３９
—
１１４０．
［１１］姜孝海，范宝春，李鸿志．基于ＡＬＥ方程的动网格
瞠口流场数值研究［Ｊ ．计算力学学报，２００８，２５
（４）：５６３—５６７．
［１２］Ｌｅｖｅｑｕｅ Ｒ Ｊ．Ｗａｖｅ ｐｒｏｐａｇａｔｉｏｎ ａｌｇｏｒｉｔｈｍｓ ｆｏｒ ｎｍｈｉ—
ｄｉｍｅｎｓｉｏｎａｌ ｈｙｐｅｒｂｏｌｉｅ ｓＹｓｔｅｎｌｓ ｆ Ｊ ．Ｊｏｕｒｎａｌ ｏｆ Ｃｏｒｎ—
ｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ，１９９７，１３１：３２７—３５３．
2024年5月26日发(作者：仇英武)
第３３卷第６期
２００９年ｌ２月
南京理工大学学报（自然科学版）
Ｊｏｕｒｎａｌ ｏｆ Ｎａｎｊｉｎｇ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｓｃｉｅｎｃｅ ａｎｄ Ｔｅｃｈｎｏｌｏｇｙ（Ｎａｔｕｒａｌ Ｓｃｉｅｎｃｅ）
Ｖ０１．３３ Ｎｏ．６
Ｄｅｃ．２００９
冲击波超压峰值的数值计算
王 杨，郭则庆，姜孝海
（南京理工大学瞬态物理国家重点实验室，江苏南京２１００９４）
摘要：为提高冲击波超压峰值的数值计算精度，该文根据二阶Ｒｏｅ格式，模拟了３种不同粗细
网格条件下的点爆炸流场。基于高精度的数值计算格式，利用冲击波的发展规律求得波阵面速
度，再由兰金一雨贡纽方程计算得到超压峰值。并与数值模拟的直接超压峰值进行了对照讨论。
结果表明随着网格加密，数值直接超压峰值不断逼近采用冲击波传播速度获得的超压峰值．而
后者几乎不变。
关键词：冲击波；Ｒｏｅ格式；冲击波阵面；兰金一雨贡纽方程；超压峰值
中图分类号：ＴＪ ０１２．２ 文章编号：１００５—９８３０（２００９）０６—０７７０—０４
Ｎｕｍｅｒｉｃａｌ Ｉｎｖｅｓｔｉｇａｔｉｏｎｓ ｏｎ Ｐｅａｋ Ｏｖｅｒｐｒｅｓｓｕｒｅ ｏｆ Ｂｌａｓｔ Ｗａｖｅ
ＷＡＮＧ Ｙａｎｇ，ＧＵＯ Ｚｅ—ｑｉｎｇ，ＪＩＡＮＧ Ｘｉａｏ—ｈａｉ
（Ｎａｔｉｏｎａｌ Ｋｅｙ Ｌａｂｏｒａｔｏｒｙ ｏｆ Ｔｒａｎｓｉｅｎｔ Ｐｈｙｓｉｃｓ，ＮＵＳＴ，Ｎａｎｊｉｎｇ ２１００９４，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ｉｎ ｏｒｄｅｒ ｔｏ ｉｍｐｒｏｖｅ ｔｈｅ ｃｏｍｐｕｔａｔｉｏｎ ｐｒｅｃｉｓｉｏｎ ｏｆ ｐｅａｋ ｏｖｅｒｐｒｅｓｓｕｒｅ ｏｆ ｔｈｅ ｂｌａｓｔ ｗａｖｅ，
ｔｈｒｅｅ ｐｏｉｎｔ ｅｘｐｌｏｓｉｏｎ ｃａｓｅｓ ａｒｅ ｓｉｍｕｌａｔｅｄ ｎｕｍｅｒｉｃａｌｌｙ ｗｉｔｈ ｄｉｆｆｅｒｅｎｔ ｆｉｎｅ ｇｒｉｄｓ ｂｙ ｕｓｉｎｇ ｔｈｅ ｓｅｃｏｎｄ－－ｏｒ－・
ｄｅｒ Ｒｏｅ ｓｃｈｅｍｅ．Ｂａｓｅｄ ｏｎ ｔｈｅ ｎｕｍｅｒｉｃａｌ ｃｏｍｐｕｔａｔｉｏｎａｌ ｓｃｈｅｍｅ ｏｆ ｔｈｅ ｈｉｇｈ—ｐｒｅｃｉｓｉｏｎ，ｔｈｅ ｂｌａｓｔ
ｐｒｏｐａｇａｔｉ’ｏｎ ｓｐｅｅｄ ｉｓ ｃａｌｃｕｌａｔｅｄ ａｃｃｏｒｄｉｎｇ ｔｏ ｔｈｅ ｄｅｖｅｌｏｐｍｅｎｔ ｌａｗ ｏｆ ａ ｂｌａｓｔ ｗａｖｅ．Ｔｈｅ ｉｎｄｉｒｅｃｔ ｐｅａｋ
ｏｖｅｒｐｒｅｓｓｕｒｅ ｉｓ ｃａｌｃｕｌａｔｅｄ ｂｙ ｕｓｉｎｇ Ｒａｎｋｉｎｅ－－Ｈｕｇｏｎｉｏｔ ｅｑｕａｔｉｏｎ ａｎｄ ｃｏｍｐａｒｅｄ ｗｉｔｈ ｔｈｅ ｄｉｒｅｃｔ ｐｅａｋ Ｏ—－
ｖｅｒｐｒｅｓｓｕｒｅ ｃａｌｃｕｌａｔｅｄ ｂｙ ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ．Ｔｈｅ ｃｏｍｐａｒａｔｉｖｅ ｒｅｓｕｌｔ ｉｎｄｉｃａｔｅｓ ｔｈａｔ：ｗｉｔｈ ｔｈｅ ｉｎ。
ｃｒｅａｓｅ ｏｆ ｔｈｅ ｍｅｓｈ ｒｅｆｉｎｅｍｅｎｔ，ｔｈｅ ｄｉｒｅｃｔ ｎｕｍｅｒｉｃａｌ ｐｅａｋ ｏｖｅｒｐｒｅｓｓｕｒｅ ｇｒａｄｕａｌｌｙ ａｐｐｒｏａｃｈｅｓ ｔｈｅ ｉｎｄｉ－
ｒｅｃｔ ｐｅａｋ ｏｖｅｒｐｒｅｓｓｕｒｅ ｏｂｔａｉｎｅｄ ｂｙ ａ ｂｌａｓｔ ｐｒｏｐａｇａｔｉｏｎ ｓｐｅｅｄ，ｗｈｉｌｅ ｔｈｅ ｌａｔｔｅｒ ｉｓ ａｌｍｏｓｔ ｃｏｎｓｔａｎｔ．
Ｋｅｙ ｗｏｒｄｓ：ｂｌａｓｔ ｗａｖｅｓ；Ｒｏｅ ｓｃｈｅｍｅ；ｓｈｏｃｋ ｆｒｏｎｔｓ；Ｒａｎｋｉｎｅ－Ｈｕｇｏｎｉｏｔ ｅｑｕａｔｉｏｎ；ｐｅａｋ ｏｖｅｒｐｒｅｓ‘
Ｓｌ１ｒＰ
火炮发射时，高温、高压火药燃气在弹丸飞出
冲击波超压峰值的测量精度与计算精度就成了本
后，自膛Ｅｌ急剧膨胀，形成膛口冲击波，对人员、装
备及周围设施形成了严重的危害。为此，各国军方
都制定了相应的安全和防护标准¨＇２ Ｊ，作为武器定
领域一直在探讨的理论与技术课题 。由于冲
击波阵面是一个间断面，因此超压峰值的测量与计
算精度受到测量系统响应与计算格式及网格的影
型与装备使用的依据之一。而冲击波的超压峰值
与冲量是该防护标准的２个基本参量。于是提高
收稿日期：２００９—０８—１０ 修回日期：２００９—１０—２７
响。近年来，一直在寻求提高精度的方法。
在试验测量方面，因冲击波超压测试系统带
基金项目：国防科技重点实验室基金
作者简介：王杨（１９７２一），女，博士生，主要研究方向：偏微分方程数值解，兵器发射理论与技术，Ｅ－ｍａｉｌ：ｗａｎｇ），一
ａｎｇｘｕｐｕ＠１６３．ｃｏｎ。
总第１６９期 王杨郭则庆姜孝海 冲击波超压峰值的数值计算 ７７１
宽有限，使得由信号响应引起的测量误差为５％
～
１０％，而这种误差一般与冲击波强度直接相关，
因此，系统实时修正及现场标定就成了保证冲击
波超压峰值测量精度的基本途径？在冲击波测量
系统标定方法的研究中，人们提出了不同的标定
方法：文献［６—８］针对不同类型的传感器提出
了其实验室标定方法和现场标定方法：美国陆军
试验规程 巾规定了利用冲击波速度标定超压
峰值的方法
对于满足物理方程的简单模型（例如不考虑
外界反应的强爆炸冲击波模型），可以得到解析
解，其超压峰值的计算是精确的 但是，描述复杂
流场冲击波（如考虑反压的弱冲击波和反射、绕
射冲击波等）运动规律的控制方程是非线性偏微
分方程组，一般多采用数值计算进行离散求解：
只有在极少情况下，事先预知激波位置，可以用
“激波装配法”得到准确的超压峰值。多数情况
下均采用“激波捕捉法”，其数值解的精度与所选
取的格式及网格粗细有关：而计算格式中加入的
人 １ 粘性，在数值求解中抹平了间断 ：选取间
断面附近不同的格点，将会得到不同的超压峰值
（文中将这种由计算格点直接获取的超压峰值称
为直接超压峰值，而由数值结果通过某种方法间
接获得的超压峰值，称为问接超压峰值） 在一
定意义上，随着网格间距的减小，直接超压峰值的
精度会有所提高。然而，随着冲击波的传播（尤
其是冲击波的远场传播），其计算区域迅速扩大，
计算量剧增，对计算机的计算速度、存储量等提出
了更高的要求；虽然存在局部加密方法（如自适
应网格等），但计算量与计算时问的矛盾仍比较
突出：如何采用较少的网格，计算 足够精度的
超压峰值就成为迫切需要解决的问题：
本文基于高精度的数值计算格式，利用冲击
波的发展规律求得波阵面速度，再由兰金一雨贡纽
方程计算得到超压峰值。根据不同网格疏密条件
下的计算结果，讨论了直接超压峰值和间接超压
峰值的变化特点＝
基本方程和计算格式
在不考虑粘性和化学反应的条件下，理想气
体满足
杀 Ｑｄ
¨ｔ）
Ｖ＋
ｆｔ）
ｄＳ＝０ （１）
式中：
Ｑ：［Ｐ ｐＵ Ｅｒ
Ｆ＝Ｑ ＋Ｇ
Ｇ＝［０ Ｐｎ ＰＵ］
ｐ、Ｕ、Ｅ分别表示流体密度、速度和单位体积的内
能，Ｓ（ｔ）表示ｔ时刻控制体Ｖ（ｔ）的表面积，ｎ表示
其外法线方向。为使方程闭合，引入状态方程
Ｐ＝ＪＤ（ —１）ｅ
式中：比热比ｙ＝１．４，ｅ为比内能。
将方程（１）离散，得
Ａ ｔ
＋
主 Ｓ：
＿一
‘‘
ｉ＝
Ｉ
０
即
＝
（ 古 Ｎ ＦｉＳ Ａ ） （２）
式中：上标ｎ表示前一时刻，ｎ＋１表示所求的时
刻。Ⅳ表示所求控制体所包含的外表面总数，Ｓ
表示第ｉ个面的面积，Ｆ 表示Ｆ在该面外法线方
向的分量，即界面流量；将控制体界面视为一维
黎曼问题，界面流量根据Ｒｏｅ方法求解
Ｆ ＝Ｆ （ＱＲ，ＱＩ）Ａ＝ＯＦ／ＯＱ
Ｆ ＝［Ｆ（Ｑ ）＋Ｆ（Ｑ ）一∑ ＩＡ，Ｉ ｅ，］／２
式中：Ｑ 和Ｑ 分别表示界面左右两侧格点值，Ａ
为Ｅｕｌｅｒ方程的Ｊａｃｏｂｉａｎ矩阵， ，Ａ ，ｅｊ分别表示
波强度矩阵 、特征向量矩阵 和右特征矩阵日
的分量。
为提高精度，引入二阶修正项
：
∑ （１一 Ａｔ ）ｗ ／２
式中：
Ｑ 一Ｑ ＝∑Ｗｉ ∑ｏｌｉ
因此，界面流量为
Ｆ ＝［Ｆ（Ｑ ）＋Ｆ（Ｑ ）一∑ ＩＡ ｌｅｊ］／２＋
∑ （１一 Ａｔ ）Ｗｉ／２
在跨声速或超声速计算中，当流场产生间断
时（如激波），其周围会产生所谓的伪震荡（Ｓｐｕｒｉ—
ＯｔＩＳ ｏｓｃｉｌｌａｔｉｏｎ） ，出现非物理解。因而需要对
其进行修正，即
Ｗｉ＝ Ｉｒ
式中：
＝（ｂ（０ ）
Ａ 如
：
Ｉ：
Ｏｔｉ Ｉ ＋ｌ Ａ ≤０
７７２ 南京理工大学学报（自然科学版） 第３３卷第６期
关于 （０）的取值，本文采用Ｍｏｎｏｔｏｎｉｚｅｄ Ｃｅｎｔｅｒｅｄ
方法，即
（０）＝ｍａｘ（０，ｒａｉｎ（（１＋０）／２，２，２０））
如前文所述，称为间接超压峰值。为了计算超压
峰值，只需在一般流场计算完毕后，再调用计算超
压峰值的子程序进行处理即可。
求得的界面流量具有二阶精度，代入式（２），即可
求得新时刻的值。对于数值解的时间步采用二阶
精度的Ｒｕｎｇｅ．Ｋｕｔｔａ法。
本文考察的点爆炸流场初始条件为：半径为
３算例与讨论
２冲击波超压峰值计算方法
根据前述数值方法获得的流场数值解，利用
波阵面上压强梯度变化最大的特征来确定波阵面
的位置，再由波阵面的位移△ｒ和相应时间间隔
△ ，可求得在此时间间隔内的平均速度。由于
一
０．１ ｍ的静止气球中充满高温、高压的高能气体，
其平均压力为２５６ｐ。，平均密度为８２ｐ。，外部为标
准条件下的大气环境，即Ｐ ＝０．１ ＭＰａ，密度为
Ｐ ＝１．１８ ｋｇ／ｉｎ 。整个计算域为半径１ １ＴＩ的球
体，外边界为外推边界。
分别模拟了网格间距为１／５００ ｍ、１／１ ０００ ＩＴＩ、
般很小（本文Ａｔ＝０．６８２５×１０ Ｓ），因此波阵
１／１ ５００ ｍ的３种情形，对应算例为Ｔ１、 、Ｔ３。
在获得数值解之后，再调用前述的计算超压峰值
的子程序分别计算间接超压峰值Ｐ （见表１）。而
直接超压峰值的取值格点如图１所示，即直接超
压峰值参考Ｐ ，Ｐ 两格点处的超压值。
面速度Ｄ可以表示为
１ ３： ：
ｔ— ●ｄ ｔ ≈——＝—Ａｔ
，一ｔ
二 （３）
、
【ｊ】
式中：ｒ。，ｒ 分别为波阵面在ｔ ，ｔ 时刻的位置。
利用激波关系式
丝：
ｐ１ ｙ＋ｌ
一
ｙ＋ｌ
式中：ｐ 、Ｍａ。＝Ｄ／ｃ 、ｃ。分别为波前气体压力、波阵
面马赫数、声速，Ｐ，为波后气体压力，进一步整理得
＝
（Ｍａ２￣－１）△ｐ＝Ｐ —Ｐ （４）
通过式（４）可以求出冲击波的超压峰值，这种由
波阵面速度和激波关系式间接获取的超压峰值，
图１ 直接超压峰值的选取示意图
表１ 算例Ｔｌ～１、３间接超压峰值与直接超压峰值列表
由表１可以看出，随着网格问距的减小，直接
为了考察直接超压峰值的精度与网格粗细的
关系，计算了算例Ｔ１、Ｔ２、Ｔ３中的直接超压峰值
Ｐ 与算例１、３的间接超压峰值Ｐ 的相对误差Ｅ 、
Ｅ，、Ｅ （见表２），随着网格的加密，相对误差越来
越小，说明直接超压峰值的精度受网格粗细的影
超压峰值Ｐ 、Ｐ 越来越逼近间接超压峰值Ｐ ，而
间接超压峰值并未受网格间距的影响，其值Ｐ 变
化很小，趋于恒定。也就是说在网格较粗的条件
下，通过间接方法得到的超压峰值就能获得较高
的精度（密网格条件下的直接超压峰值）。
响较大。为了得到较精确的直接超压峰值，必须
总第１６９期 王 杨 郭则庆 姜孝海 冲击波超压峰值的数值计算 ７７３
加密网格。但必须注意到，由于加密网格而引起
的计算量、存储量以及计算时间等将会大大增加：
可见，在本文的研究范围内，采用间接方法计算超
压峰值较直接方法（加密网格）更有现实意义和
可行性
表２算例Ｔ１～１、３直接超压峰值与 间接
超压峰值的相对误差列表
图２为算例１、３的在不同时刻的压力和径向
速发典型分布图：从图２（ａ）可见，初始冲击波形
成的同时，膨胀波从接触面开始朝球心内部运动
球心内部点的压力初始保持定值，膨胀波经过后，
压力降低，形成第二道激波，见图２（ｂ）：图２（ｃ）
展示的是典型的压力图，主冲击波后是负相区：
１・４
１ ２
日
皇０．８
运
Ｏ４
０
１５０ ２ｏ０ ２５０ ３Ｏ０ ３５０
ｍｍ
（ａ）ｎ＝１
ｆｂ）”，＝６
３５
２
３０
２５
１
２０ ０
１５
一
１
１０
－
０５
２
６００ ８００
ｒ／ｎｌｒｆｌ
（ｃ）ｎ ＝１６
图２不同时刻的压力与速度分布曲线
（ｎ，表示时间步数）
４结束语
本文采用二阶精度Ｒｏｅ格式，模拟了３种不同
粗细网格的点爆炸流场：利用激波关系式和波阵
面速度间接获得了冲击波的超压峰值 与直接超
压峰值的对比分析表明：采用间接计算方法，可以
在粗网格条件下获得比较准确的初始冲击波超压
峰值 这在一定意义上解决了计算区域大、网格数
量多造成的计算机时长与精度之间的矛盾：
参考文献
『１］ ＭＩＩ 一ＳＴＤ １４７４（ＭＩ）一１９７３，Ｎｏｉｓｅ Ｌｉｍｉｔｓ ｆｏｒ ａｒｍ ｙ ｍａ—
ｔｅｒｉｅｌ［ｓ］．
『２］ ＧＪＢ １１５８—９１，炮口冲击波对人员非听觉器官损伤
的安全限值［ｓ］．
『３ １ Ｋｉｎｎｅｖ Ｇ Ｆ．Ｇｒａｈａｍ Ｋ ｊ．Ｅｘｐｌｏｓｉｖｅ ｓｈｏｃｋｓ ｉｎ ａｉｒ
［Ｍ］．Ｂｅｒｌｉｎ ａｎｄ Ｎｅｗ Ｙｏｒｋ：Ｓｐｒｉｎｇｅｒ—Ｖｅｒｌａｇ，１ ９８５．
［４］ Ｖａｎｄｅｒｓｔｒａｅｔｅｎ Ｂ，ｌ￣ｅｆｅｂｖｒｅ Ｍ，Ｂｅｒｇｈｍａｎｓ Ｊ．Ａ ｓｉｍ‘
ｐｉｅ ｂｌａｓｔ ｗａｖｅ ｍｏｄｅｌ ｆｏｒ ｂｕｒｓｔｉｎｇ ｓｐｈｅｒｅｓ ｂａｓｅｄ Ｏｎ ｎＵ—
ｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ『Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ｈａｚａｒｄｏｕｓ Ｍａｔｅｒｉ—
ａｌｓ，１９９６，４６（２）：１４５—１５７．
［５］李鸿志．炮口制退器科研阶段报告［Ｒ 南京：华东
工程学院，１９７８．１０．
［６］ 张远平，池家春．爆炸冲击波压力传感器灵敏度的
动态标定及测试技术研究［Ａ］．第四届全 爆炸力
学实验技术学术会议论文集［Ｃ ．绵阳：中国力学
学会，２００６．３４９—３５３．
［７］ 崔海涛，刘庆明．冲击波压力传感器测试系统的动
态标定［Ｊ］．流体力学实验与测量，２００４，１８（１）：９２
—
９６．
［８］ 王等旺，张德志，李焰，等．冲击波压扦式传感器测
试系统的动态标定［Ａ］．第十二届全国激波与激
波管学术会议论文集［Ｃ ．洛阳：中国力学学会，
２００６．３６３—３６６．
［９］ ＭＴＰ ４—３—８２２，美陆军试验鉴定部兵器试验规程一
冲击波压力测量（电子的）技术试验［Ｓ ．
ｆ ｌ０ １ Ｍｕｒｄｏｃｋ Ｊ Ｗ．Ｓｈｏｃｋ—ｗａｖｅ ｉｎｔｅｒａｃｔｉｏｎ Ｍｔｈ ｔｗｏ．ｄｉｍｅｎ—
ｓｉｏｎａｌ ｂｏｄｉｅｓ［Ｊ］．ＡＩＡＡ Ｊｏｕｒｎａｌ，１９７５．１３（９），１１３９
—
１１４０．
［１１］姜孝海，范宝春，李鸿志．基于ＡＬＥ方程的动网格
瞠口流场数值研究［Ｊ ．计算力学学报，２００８，２５
（４）：５６３—５６７．
［１２］Ｌｅｖｅｑｕｅ Ｒ Ｊ．Ｗａｖｅ ｐｒｏｐａｇａｔｉｏｎ ａｌｇｏｒｉｔｈｍｓ ｆｏｒ ｎｍｈｉ—
ｄｉｍｅｎｓｉｏｎａｌ ｈｙｐｅｒｂｏｌｉｅ ｓＹｓｔｅｎｌｓ ｆ Ｊ ．Ｊｏｕｒｎａｌ ｏｆ Ｃｏｒｎ—
ｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ，１９９７，１３１：３２７—３５３．