《凝聚态炸药爆轰数值模拟》由会员分享，可在线阅读，更多相关《凝聚态炸药爆轰数值模拟（5页珍藏版）》请在金锄头文库上搜索。
1、凝聚态炸药爆轰数值模拟程俊霞( 北京应用物理与计算数学研究所北京)蔫要：本文以工程数值模拟方法为主，回顾了凝聚态炸药数值模拟的模型，反应率和状态方程。关键词：凝聚态炸药反应率状态方程1 引言1 9 世纪8 0 年代初B e r t h e l o t 和V i e i l l e 及M a l l a r d 和L e C h a t e l i e r 在观察管道中燃烧火焰的传播过程中发现了爆轰波。如今一百多年过去了，虽然炸药广泛应用，但人们对爆轰的理解仍然非常有限。爆轰的理论模型可分为理想和非理想爆轰模型。c J 模型和Z N D 模型是理想一维爆轰模型。C h a p m a n 和J
2、o u g u e t 分别于1 8 9 9 年和1 9 0 4 年各自独立提出了理想爆轰模型，后来称之为c J 模型。c J 模型是一维模型，驱动爆轰的冲击波是平面的，并且不考虑反应区，认为反应是瞬间完成的，反应完成的状态满足C J 条件。理想爆轰速度定义为与化学释放能和状态方程的函数，但是与能量释放速度无关。二战期间，Z e l d o v i c h ，V o n N e u m a n n 和D o r i n g 提出了带化学反应区的一维稳态爆轰模型，称之为Z N D 模 型。在冲击波压缩后开始反应，完全反应的平衡状态对应于c J 理论的c J 点。因此，也被称为理想爆轰模型。当药柱
3、直径有限时，爆轰波阵面总是弯曲的，反应产物飞散，在爆轰反应区后反应并不完全，边界处的爆速变慢，当药柱的直径小到一定程度时，爆轰熄灭称此药柱直径为临界直径。常常称这种模型为非理想爆轰模型弘J 。2 爆轰模型爆轰波是冲击波和化学反应的耦合作用。在纳秒时间内，在很薄的反应区内，物质温度升高3 0 0 0 k 5 0 0 0 k ，压力升至几十万大气压。虽然理论和试验方面都朝着揭开这个极快而剧烈的神秘现象面纱而努力，但是直到目前为止，我们对此仍知之甚少。随着凝聚态炸药日益广泛的应用，尤其是在武器研究方面，爆轰的数值模拟显得颇为重要。但是，由于爆轰过程的时间和空间的巨大跨度，使爆轰的数值模拟非常困难。爆
4、轰的反应区一般在0 0 1 1 0毫米【l 】，而我们所研究的问题是在1 0 0 毫米的量级。高能炸药的反应区中化学反应时间一般在皮秒量级，而我们计算的问题爆轰过程也在几十微秒这么长的时间。鉴于爆轰过程的微观特性，有人采用分子动力学方法1 3 , 4 模拟爆轰过程，研究爆轰过程的细节，试图得到爆轰反应机制的信息，但是这种方法还处于起步阶段，化学反应模型还停留在简单反应模型阶段。另外，分子动力学的计算量太大，很难在工程中实际应用。爆轰的分子动力学模拟对我们进一步认识爆轰，发展爆轰理论具有参考意义，但与实际工程应用的距离较远。工程上要求爆轰数值模拟求速度快，但叉能反映非理想爆轰的基本特征。基于这个
5、要求提出了一些爆轰模型。在工程设计中，经典模型由于其求解速度快并且结果基本能与试验符合，因此一直应用至今。但是，提高数值模拟的精度的追求和要求一直没有改变。由于一维理想爆轰在解释爆轰熄灭、边界效应等现象时显得无能为力，早在二次世界大战期间，人们开始注意到二维爆轰。真空边界处压力和爆速的降低，汇聚爆轰的高爆速等效应都是反应区作用的结果。在一般工程应用中，特别是当我们计算爆轰波传播问题中，常常忽略反应区的影响。当反应区的宽度比所研究的问题特征尺度小的多时，反应区效应不重要，忽略反应区带来的误差不大。但是对于反应区厚度约毫米量级的钝感炸药来说，忽略反应区将不是一个可靠的近似。此时，能刻画反应区效应的
6、工程模拟方法亟待出现。1 6非理想爆轰的研究开始的很早。D e v o n s h i r e B 3 1 研究了非理想爆轰，并且推出了著名的粒子速度梯度方程，正确地推断了当反应粒子的相对速度为当地声速时，爆轰将会终止。W o o d K i r k w o o d 5 1 研究了稳定二维轴对称爆轰，提出了轴摄动理论。扩展Z N D 模型描述稳定冲击波后的轴向流动。他们导出了与R H 关系相似的声速点上的能量释放率和流发散的积分守恒方程。根据W o o d K i r k w o o d ，稳定爆轰速度随波阵面曲率的增加而减小。W o o d K i r k w o o d 理论并不限定产物反
7、应物的状态方程和反应率的形式。给定两相状态方程和反应率，对称轴处的曲率给定的条件下，W o o d - K i r k w o o d 理论就给出了爆轰速度，反应区宽度，压力分布，温度，粒子速度和反应度。这个理论扩展到多维情况可能具有很大的困难，直到1 9 8 1 年B d z i l 9 1 才提出了二维理论。B d z i l 和S t e w a r t 在B d z i l 的基础上提出了渐进方法D S D 方法【l0 1 。他们在爆轰波反应区宽度远小于爆轰波曲率半径的极限下，用简单的方法确定爆轰波的发展方程和反应区中的流。反应欧拉方程的直接数值模拟确认了D S D 理论的主要结论1
8、1 1 j ，同时也揭示了二维反应区结构的其它特征。他们考虑了凸冲击波后的反应区的径向结构。对于无约束边界和弱约束边界，冲击波和声速面相交于炸药柱的侧边界上。但是D S D 理论的基础是冲击波的曲率半径远大于药柱的半径，且反应区的宽度远小于药柱的半径。边界条件采用半经验参数。B d z i l 分析了小弯曲冲击波阵面和薄反应区的情况，得到了描述爆轰速度和药柱直径之间关系的冲击波形状的微分方程。B d z i l 和S t e w a r t 在此基础上引入了时间项，提出了一个半经验的理论，称之为D S D 方法，主要用以研究爆轰波在复杂几何形状中的传播过程。从其本质而言，D S D 方法不适合
9、模拟爆轰熄灭问题。D S D 方法是基于爆轰波的法向速度见与爆轰波的曲率r 之间的关系，它模拟的主要是曲率效应对爆轰发展过程的影响。爆轰波的速度和曲率存在一定的关系。K e n n e d y 1 2 1 也分析了三种不同炸药非理想爆轰的一系列的试验数据，结果显示，在最不理想的爆轰中，也存在见，K 的关系。其它两种炸药支持了Y a o & D S S t e w a r t ( 1 9 9 6 ) 包括爆轰加速度见的乜，彭关系。目前D S D 方法已经应用于工程设计中。3 反应率人们对炸药的反应率、反应物产物的状态方程也进行了广泛而长期的研究。由于炸药化学反应复杂而迅速，因此还处于唯象模型的阶
10、段，反应率只是一些热力学量的函数，需要通过试验进行参数的标定。常用的定常的反应率有W i l k i n s 函数和C J 比体积燃烧函数，因为程序计算时需要给出c J 模型的伪反应区而给出的。非定常爆轰的反应率函数有：A r h e n i u s 反应率、双A r h e n i u s 反应率、F o r e s tF i r e燃烧函数、C o c h r a n 起爆函数、L e e - - T a r v e r 点火成长模型和T a n g 统一模型速率方程。W i l k i n s 函数【8 】是使用最早的一种反应率函数，它与爆速和网格大小有关。A =+ 筹( 3 1 )其中
11、欲为网格的特征长度，一般取= 1 ，口 l ，Z 为根据网格与起爆点的位置根据惠更斯原理得到的起爆时间。C J 比体积燃烧函数定义为1 7R D竺加+II 一 一0qIq墨欲f D口五= f 旦兰l( 3 - 2 ) 一吁m 为可调参数。由于其形式简单，实际使用较多。以上两种反应率用于定常爆轰的模拟。下面介绍在非定常爆轰模拟中经常采用的反应率。A r r h e n i u s 反应率也是一种唯象的反应率，但是与欧拉方程相耦合，常常给人以准确模拟的错觉。在B d z i l 的D S D 方法确认中，称欧拉方程和A r r h e n i u s 反应率的求解为直接数值模拟。它的形式为i d
12、2 = 2 Z e l 冲( _ 寺B 3 ，A U II J J ， 破1 lR 7J、其中Z ，E ，R 均为常数，T 为温度。A r r h e n i u s 反应率与定常爆轰反应率最大的不同除了与时间相关外，还与温度相关。与此相似的是双A r r h e n i u s 反应率，也成为两步模型。警一毛p 唧( 一鲁)警吐p z 徊( 一鲁) 卜五，唧( 一警)。4 F o r e s tf i r e 反应率为压力的函数，是以试验测出的炸药物质带化学反应的雨贡纽关系拟合得到的函数譬：( N e x Da )( 3 35 一)= 一九I7；口I( )者、智“7、。它最大的缺点是各项没有
13、明确的物理意义，只是一个拟合的关系式，也或许是这样的物理意义还没有被发现。C o c h r a n 起爆函数它的每一项具有一定的物理意义，并反映一定的爆轰起爆和发展的机制，它的形式为i d 2 = ( 1 一旯) q p ”+ 卸1( 3 6 )蚴，吨，m ，n 为试验确定的常数。L e e T a r v e r 反应率函数，也成为点火成长模型，在非定常爆轰中较为常用，其形式为皇= I ( 1 一兄) 。( V o v 1 ) 7 + G ( 1 一A ) 。允y p 2( 3 7 )第一项为点火项，第二项为成长项。T a n g 统一模型速率模型是将炸药的起爆、爆轰机制分解：热点内炸药的
14、分解，热点周围炸药的预热与中间产物的形成，炸药化学反应分解为快慢反应。4 状态方程对于爆轰问题，产物和反应物的状态方程需要分别给出，反应区内的状态方程采用二者的组合。与产物的状态方程相比，反应物的状态方程研究较少。当计算中需要反应物的状态方程时，一般用G r u n e i s e n 状态方程，或者J W L 状态方程。直到目前，大多状态方程限制在求解无粘，无热传导，无反应和无辐射的流体力学方程中，因此，不完全状态方程已经足够了。不完全状态方程不包含温度和熵，对于求解理想炸药作为驱动的问题足够了。计算非理想炸药的反应区，D D T ，临界直径，部分反应等1 8问题时需要求解完全状态方程给出与
15、反应率所依赖的温度。在产物的状态方程中，最简单的是多方气体状态方程。 P = R T V( 4 1 ) 其中瞒比体积，R 为气体常数，功温度。由于反应产物的高温高压高密度，理想气体状态方程是不适用的。第一个简单的描述产物的经验状态方程是B K w 状态方程( B e c k e r ，1 9 2 1 ；K i s t i a k o w s k y & W i l s o n ，1 9 4 1 ) 。在经验状态方程中，应用的最成功的是J W L 状态方程( J o n e s W i l k i n s L e e ，1 9 6 8 ) ，至今仍然有广泛的应用。J W L 状态方程描述了产物的
16、绝热膨胀过程，其参数由圆筒试验给定。p 刊( 卜喝即一护喝等其中瞒相对体积，V = P o p 。在实际应用中，反应物的状态方程与产物的状态方程类似，也采用J W L 状态方程形式。混合物采用二者的组合。国为G n m e i s e n 系数，认为是一个常数。P M V 状态方程q 是描述反应物的完全状态方程，S L 茭j H e l m h o l t z 由能4 口，y ) 的形式，由压力关系导出尸= ( p l p o ) 2 P 2 + ( p l p o ) 3 伤+ ( p l p o ) 4 肋+ O T( 4 3 )定容比热为G = Q + C l ( 瓦丁) + C 2 ( T o 丁) 2( 4 4 )其中p ，C 。，0 均为常数。反应物的孙状态方程需正确地与试验冲击雨贡纽线相对应。W i l l i a m b u r g ( w B G ) 状态方程是描述产物的完全状态方程，以比体积和熵为自变量。 E = 丽P V( 4 5
