Monte-Carlo analysis of the projectile fragments from cylindrical tank boiling liquid expanding vapor explosion accident
-
摘要: 在动力学分析的基础上,将二维体系内的碎片轨迹方程扩展到三维体系内,将风这一影响因素考 虑在内。以墨西哥事故中水平圆柱储罐发生沸腾液体扩展为蒸气云爆炸(BLEVE)为例,利用Monte-Carlo法 模拟出碎片的轨迹曲线,计算得到有风和无风情况下碎片抛射距离的概率分布,得到一般情况下风对碎片危 害的影响较小。计算了无风情况下碎片碰撞目标容器的概率,得到碰撞概率随罐间距的增大而呈现严格的指 数衰减趋势,通过此关系式可以计算出符合安全要求的罐间距。研究结果对于提高储罐的安全性、缓解和控 制碎片产生的风险具有一定的指导意义。
-
关键词:
- 爆炸力学 /
- 沸腾液体扩展为蒸气云爆炸 /
- Monte-Carlo法 /
- 碎片 /
- 柱形罐
Abstract: Thegenerationandflightoffragmentsarestochasticandcancausemoreseverehazards. Studyingonthedominoeffectstudycausedbythefragmentsfromexplosionishotspotanddifficulty. Basedondynamicanalysis,thepresentworkextendedthetwo-dimensionalfragmenttrajectoryequationintothree- dimensionalsystem,andtookwindintoaccount.ThepresentanalysistookacylindricaltankinMexicoCitywithboilingliquidexpandingvaporexplosionasanexample, andthetrajectoryoffragmentwassimulatedbyMonte- Carlomethod.Seenfromtheprobabilitiesoffragmentranges tobereachedwith/withoutwind,windhaslessimpactonthehazardcausedbyfragments.Theimpactprobabilityofafragmentonatargetwithnowindwasattained, impactprobabilitywiththetank showedtheincreaseofdistancebetweenthestrictexponentialdecaytrend,andthroughthisrelationship, thetankspacingincompliancewithsafetyrequirementscanbecalculated.Thepaperprovidesa methodwhichisofgreatsignificancefortheriskassessmentandcontroloffragmentsfromtheexplosionofcylindricaltanks. -
电磁轨道炮是一种利用高功率脉冲电流产生的巨大电磁推力把弹丸加速到超高速的动能武器,在军事上有着广泛应用前景。电磁轨道炮走向应用必须要攻克发射器的使用寿命这项关键技术,延长发射器寿命则需要解决电弧烧蚀和刨削等问题[1]。刨削是材料在高速滑动接触下出现的一种损伤现象,在火箭橇、轻气炮及电磁轨道炮的实验中均有出现[2-3]。轨道炮中的刨削主要带来以下危害:降低能量利用率,缩短轨道寿命,影响射弹精度等。轨道炮实验中观察到的典型刨坑形貌多为泪滴状,如图 1所示[4],刨坑尖部朝炮尾、弧部指向炮口。刨坑的生长,从一个起始点开始,沿电枢的运动方向逐渐变宽变深,坑内表面凹凸不平,呈现熔融的表面状态,刨坑最终收敛于弧形的尾边,尾边上常有熔融物堆积成的凸缘。C.Persad[5]在轨道炮实验中测量出发生刨削位置的轨道背面表面高度差达到50 μm,而未发生刨削位置的高度差不超过5 μm,说明刨削的产生过程伴随着很强的冲击。T.J.Watt[6]通过实验研究对比了平整轨道和表面带凹缝的轨道对刨削临界速度的影响,发现轨道表面宏观不平度对产生刨削的临界速度影响不大,而对刨坑的形貌影响较大。F.Stefani[4]对不同材料的刨削临界速度进行了实验,认为枢轨中较硬材料的维氏硬度与刨削临界速度下的冲击压强近似满足线性关系。
电枢与轨道超高速滑动接触界面上刨削现象的形成是一个瞬态过程,持续时间通常只有几十微秒,加之该过程处于极端工况下,以目前的实验水平对轨道炮刨削的产生过程进行实时测量基本不可能。因此,若能对轨道刨削现象的演变过程进行合理的复现,刨削的形成机理及相关规律的研究便可迎刃而解。本文中基于轨道炮结构特点及冲击热力学理论,通过物质点法以三维数值模拟的方式研究轨道刨削的产生机理,展现刨削的演化与形成过程,分析刨削产生的极端物理环境,并在此基础上,进一步计算和分析颗粒大小、冲击速度等因素对刨削产生的影响。
1. 轨道刨削机理及数值研究方法分析
J.P.Barber[7]首次分析了轨道炮的刨削产生机理,参照前人在火箭橇中的研究成果,把轨道炮刨削的产生归因于轨道微颗粒的冲击作用,认为刨削是在超高速滑动接触时碰撞产生的冲击应力超过材料的极限强度后产生的。T.J.Watt[6]的实验也表明,磨损以及刨削,均起因于轨道表面的缺陷,这些缺陷通常是机械加工或打磨时留下的痕迹。L.M.Barker[8]首次将轨道炮刨削机理用平行冲击热力学理论加以解释,电枢相对轨道平行滑动,由于轨道表面的不平度或细小颗粒,引起电枢与轨道间的微小冲击,并超过材料的应力屈服强度,枢轨材料发生塑性变形,最终在轨道上形成刨坑。实际上,轨道表面微颗粒除机械加工因素外,还可来自以下方面:电枢运动过程中发生折断等破坏,产生的碎片或残渣附着于轨道表面;电枢局部积聚热量过高,使得电枢材料熔化掉落于轨道表面,凝固成细小颗粒[9];由于风沙等环境因素,也可导致轨道表面附有细小沙粒。
数值模拟方面,L.M.Barker[10]与D.L.Bourell[11]基于微颗粒诱发的刨削模型,分别用CTH和EPIC等流体动力学程序对轨道炮刨削过程进行了二维数值模拟,对刨削受材料、速度等变量的影响进行了计算,但结果中并未给出完整的刨坑形貌。金龙文等[9]、Zhu Rengui等[12]、刘峰等[13]利用LS-DYNA和ABAQUS等商业有限元软件,也对轨道炮刨削的影响因素进行了数值研究,分析得出刨坑的形貌受电枢冲击前角、颗粒尺寸、冲击速度及温度等变量影响较大,枢轨界面载荷对刨削影响较小。
目前的数值模拟多为二维,而轨道刨削是一个三维扩展问题,二维计算结果不能完整反映出刨坑形貌的扩展过程以及受计算参量变化的影响。此外,枢轨高速滑动界面上的刨削现象是局域化的,变形主要集中在接触界面附近薄层内。使用有限元法时,若想得到枢轨界面的碰撞作用和刨坑的形成过程,一是需要在界面附近划分细密的网格,势必会导致计算量增大,计算时长增加;二是超高速刨削现象涉及材料的非线性大变形破坏,有限元法计算时容易因网格畸变过大从而大幅缩短时间步长,为保证计算效率,需要进行网格重分或删除失效单元。物质点法[14-15]是粒子类无网格方法的一种,兼具拉格朗日法与欧拉法的优点。将物质区域离散为物质点,在规则的背景网格上进行动量积分,从而避免了网格畸变问题,计算时间步长也基本固定,适合于超高速刨削这类极端现象的模拟。
2. 计算模型及物质点法程序实现
2.1 计算模型及材料参数
根据上述刨削产生机理的分析,采用物质点法建立了轨道微颗粒诱发的刨削模型(图 2)。由于结构的对称性,沿电枢的对称面建立1/2计算模型。电枢对称面上施加对称约束条件,模型中略去绝缘支撑材料,改为在电枢侧面施加z向位移约束。U形电枢长30 mm,轨道宽30 mm,厚10 mm,轨道表面有一半径0.5 mm的半球形颗粒。整个模型物质域以627 792个质点进行离散,其中电枢171 280个质点,轨道454 400个质点,颗粒2 112个质点。
电枢材料为7075铝,轨道材料为无氧铜,材料本构为Johnson-Cook模型。J-C模型中考虑到应变硬化、应变率硬化及温度软化效应,流动应力表示为
σy=(A+Bεpn)(1+Cln˙εp˙ε0)[1−T∗m] (1) 式中:A是参考温度与参考应变率下的屈服强度,B、n是材料应变硬化参数;C是应变率强化参数,εp是等效塑性应变,ε0是参考应变;T*=(T-Tr)/(Tm-Tr),m是热软化参数,Tr是参考温度,Tm是材料熔化温度。
考虑损伤的破坏应变εp, f表示为
εp,f=(D1+D2eD3σ∗)(1+D4ln˙εp˙ε0)(1+D5T∗) (2) 式中:σ*=σm/σeq为应力三轴度,σm是静水应力,σeq是von Mises等效应力。D1~D5是材料断裂失效参数。当质点损伤量D=∑(Δεp/εp, f)达到1时,材料断裂失效,已失效的质点偏应力置为零,且不能受拉。
刨削的形成时间通常只有几十微秒,可视为绝热过程。绝热过程塑性变形引起的温升为
ΔT=∫εp0χσeqρcpdεp (3) 式中:ρ为密度;cp为定压比热;χ表示塑性功转化为热量的比例,一般取0.85~0.9。
由于刨削涉及到材料在冲击高压下的热力学行为,因此还须用相应的状态方程进行描述,本文采用Mie-Grüneisen状态方程。该状态方程定义材料所受压力为
p={ρ0c20μ[1+(1−γ0/2)μ][1−(s−1)μ]2+γ0ρ0eμ⩾ (4) 式中:c0是压力为零时的材料声速,s是冲击波速us和粒子速度up关系曲线的斜率系数,γ0是压力为零时的Grüneisen参数,ρ0为材料初始密度,e为比内能,体积变化率μ=ρ/ρ0-1。电枢与轨道的材料参数见表 1,其中:E为杨氏模量,ν为泊松比。
表 1 电枢与轨道材料参数Table 1. Material parameters of armature and rail材料 ρ/(kg·m-3) E/GPa ν cp/(J·kg-1·K-1) A/MPa B/MPa n C m Tm/K Tr/K {\dot \varepsilon }0/s-1 χ D1 D2 D3 D4 D5 c0/(m·s-1) s γ0 7075铝 2 810 71 0.33 960 369 684 0.73 0.008 3 1.7 933 293 1 0.9 0.13 0.13 -1.5 0.011 0 5 350 1.34 2 无氧铜 8 960 124 0.34 383 90 292 0.31 0.025 1.09 1 356 293 1 0.9 0.54 4 2 0.014 1.12 3 940 1.49 2 2.2 物质点法程序实现
冲击爆炸等问题中常使用更新拉格朗日格式,其考虑热量交换的控制方程为
质量方程\;\;\;\;\;\frac{{{\text{d}}\rho }}{{{\text{d}}t}} + \rho \frac{{\partial {v_k}}}{{\partial {x_k}}} = 0 (5) 动量方程\;\;\;\;\;\frac{{\partial {\sigma _{ij}}}}{{\partial {x_j}}} + \rho {b_i} = \rho {{\ddot u}_i} (6) 能量方程\;\;\;\;\;\rho \dot e = \rho {q_{\text{s}}} + \frac{\partial }{{\partial {x_i}}}(k\frac{{\partial T}}{{\partial {x_i}}}) + {s_{ij}}{{\dot \varepsilon }_{ij}} - (p + q){{\dot \varepsilon }_{kk}} (7) 边界条件\ \ \ \ \ ({{n}_{j}}{{\sigma }_{ij}})\left| _{{{\mathit{\Gamma} }_{\text{t}}}} \right.={{{\bar{t}}}_{i}},\ \ \ \ \dot{u}\left| _{{{\mathit{\Gamma} }_{\text{u}}}} \right.={{{\bar{\dot{u}}}}_{i}} (8) 式中:i, j, k=1, 2, 3分别代表直角坐标系中的三个方向;ρ是密度,bi是物体所受单位体积力,{\ddot u}是加速度,σij是Cauchy应力,偏应力张量sij=σij+pδij,p是压力,e是系统的比内能,qs是单位质量热源,T为温度,Γt和Γu分别是给定面力边界和位移边界,nj是边界At的外法线单位矢量,q是人工体积黏性力。
物质点法求解格式也是从微分方程的弱形式出发,取虚位移δui作为权函数,式(6)和式(8)的等效积分弱形式即为
\int\limits_\mathit{\Omega} {\rho \ddot u_i\delta {u_i}{\text{d}}V} + \int\limits_\mathit{\Omega} {\rho \sigma _{ij}^*\delta \frac{{\partial {u_i}}}{{\partial {x_j}}}{\text{d}}V - \int\limits_\mathit{\Omega} {\rho {b_i}\delta {u_i}{\text{d}}V} } - \int\limits_{{\mathit{\Gamma} _t}} {\rho \bar t_i^*\delta {u_i}{\text{d}}A} = 0 (9) 式中:\sigma _{ij}^*=σij/ρ是比应力,\bar t_i^* = {{\bar t}_i}/\rho 是比边界力,Ω为求解空间。
物质点法将连续体离散为一系列质点,因此连续体的密度可近似为
\rho ({x_i}) = \sum\limits_{w = 1}^W {{m_w}\delta ({x_i} - {x_{i, w}})} (10) 式中:W是质点总数,mw是质点w的质量,δ是Dirac Delta函数,xi, w是质点w的坐标。将式(10)代入虚功方程(9)中,可将其转化为求和的形式
\begin{gathered} \sum\limits_{w = 1}^{{n_w}} {({m_w}{{\ddot u}_{i,w}}\delta {u_{i,w}})} + \sum\limits_{w = 1}^{{n_w}} {\left( {{m_w}\sigma _{ij,w}^*\delta \frac{{\partial {u_{i,w}}}}{{\partial {x_j}}}} \right)} - \sum\limits_{w = 1}^W {({m_w}{b_{i,w}}\delta {u_{i,w}})} - \hfill \\ \sum\limits_{w = 1}^W {({m_p}\bar t_{i,w}^*{h^{ - 1}}\delta {u_{i,w}})} = 0 \hfill \\ \end{gathered} (11) 式中:下角标w表示物理量描述质点w,h是为了将式(9)等号左边最后一项边界积分转化为体积分而引入的假想边界层厚度。
物质点法中物体的信息均由各质点携带,为了在背景网格上求解tn时刻的动量方程,需要将质点在tn时刻的质量、动量、应力、体力和面力等信息映射到背景网格,利用背景网格结点上的有限元形函数NI(xi)实现质点和背景网格结点之间信息的映射,第I个网格结点在i方向的动量为
{P_{i, I}} = \sum\limits_{w = 1}^W {({m_w}{N_{I, w}}{{\dot u}_{i, I}})} (12) 结点内力和结点外力分别为
f_{i, I}^{{\text{in}}\;{\text{t}}} = - \sum\limits_{w = 1}^W {\left( {\frac{{\partial {N_{I, w}}}}{{\partial {x_j}}}{\delta _{ij, w}}\frac{{{m_w}}}{{{\rho _w}}}} \right)} (13) f_{i, I}^{{\text{ext}}} = \sum\limits_{w = 1}^W {({m_w}{N_{I, w}}{b_{i, w}})} + \sum\limits_{w = 1}^W {({N_{I, w}}{{\bar t}_{i, w}}{h^{ - 1}}\frac{{{m_w}}}{{{\rho _w}}})} (14) 式中:σij, w为质点w的应力,可以利用本构方程由质点w的应变增量和旋量增量得到。
背景网格结点的动量方程可表示为
{{\dot P}_{i, I}} = f_{i, I}^{{\text{in}}\;{\text{t}}} + f_{i, I}^{{\text{ext}}}, \;\;\;\;{x_I} \notin {\mathit{\Gamma} _{\text{u}}} (15) 在背景网格上利用时间积分求解各结点的动量方程后,再将网格结点的速度变化量和位置变化量映射回质点,以更新质点的速度和位置。
3. 结果与讨论
3.1 刨削过程模拟
以一个典型算例进行分析,给电枢赋予1.5 km/s的初始速度,与轨道上的半径为0.5 mm半球形颗粒发生平行冲击,整个作用时间约30 μs。图 3为刨坑形成过程中Mises等效应力变化云图(x为电枢运动方向,y为枢轨接触面法向,截取模型中间对称面观测,为便于查看刨坑形成过程,后处理时已破坏的质点未予显示),应力波从电枢与颗粒的碰撞点开始向周围迅速扩散,冲击核心区内应变梯度很大,材料产生大变形流动,应变率最高时可达105 s-1。电枢和轨道材料上的Mises等效应力分别达到了900和450 MPa,冲击瞬间轨道微颗粒受到的压力可至几十吉帕,远高于材料静态屈服强度,材料在急剧压缩下表现近似流体状态。冲击塑性功产生的热量来不及扩散,使得电枢和轨道局部温度升高了300 K以上,进一步加剧了冲击核心区枢轨材料的热软化及流动变形。图 4是冲击过程中枢轨接触面上的作用力随时间的变化曲线,在电枢与微颗粒发生冲击作用的时间内,接触力发生跃变。x及y方向力较大,峰值可达15 kN。y方向的力使得轨道表面材料被急剧压缩,刨坑内部的材料密度升高,即轨道材料被硬化。文献[4]对轨道刨坑区域的材料硬度测量也显示,刨坑内材料硬度增加明显,表明刨削发生时在枢轨接触面法向的轨道材料塑性压缩变形效应显著。
图 5是轨道微颗粒内选取的2个质点w1、w2在冲击发生后的速度(vx、vy)变化曲线,质点在x和y方向的动量均迅速增长(负号表示y轴反向),说明轨道表面被斜角度侵彻,引起该处材料的流动与扩散,这2个质点在17 μs后均发生了破坏。图 6是冲击结束后,轨道表面形成的泪滴状刨坑,尖端指向轨道尾部,弧端指向炮口,刨坑长12 mm,宽8 mm,深1.2 mm。计算模拟的结果与轨道炮实验中的泪滴状刨坑形貌(图 1)符合较好。与二维计算[10-12]相比,三维模拟能比较清晰地反映出材料塑性波从冲击核心区域向空间三维的扩展过程及最终形成的泪滴状刨坑形貌。
对刨削形成过程的数值模拟表明,轨道刨削是在一个瞬态的局域高压扰动下产生的。轨道微颗粒诱发的瞬态局域冲击使得枢轨界面材料发生高应变率的塑性流动和强化等热力学行为。从上游而来的电枢与轨道表面的局域高速冲击产生了瞬时的能量交换,使得受冲击的轨道材料变成了高热高压的金属流,并以斜角度侵彻下游的轨道表面,轨道材料发生塑性流动。塑性变形从冲击扰动点开始,向下游沿3个方向逐渐扩展,最终形成了泪滴状的刨坑。
3.2 影响因素分析
保持枢轨结构及材料不变,电枢仍赋予1.5 km/s的初始速度,逐渐减小微颗粒的半径尺寸,分析其对刨削的影响。结果显示,随着颗粒尺寸的减小,冲击形成的刨坑尺寸也逐渐减小,颗粒半径为0.2 mm时形成的泪滴状刨坑如图 6(b)所示。计算结果随颗粒半径的进一步减小呈现出收敛特征,当颗粒半径减小到0.1 mm以下后,轨道表面泪滴状刨坑逐渐消失,相应的,电枢的损伤破坏也逐渐减轻。
改变电枢与轨道颗粒发生冲击的速度,分析轨道刨削受冲击速度的影响。结果显示,随着冲击速度的提高,形成的刨坑变深变宽,电枢的损伤面积也变大,如图 7所示。计算得到的7075铝电枢与无氧铜轨道发生刨削的临界速度约为1.4 km/s,冲击速度低于1.4 km/s后,轨道上不再形成泪滴状刨坑,而是出现了磨损擦伤,如图 8所示,与刨坑不同,磨损形成的沟槽较浅,且沿生长方向基本保持恒定的宽度和深度。随着电枢冲击速度的减小,磨损也逐渐减轻,直至消失。这是由于刨削的形成是轨道表面受电枢高速冲击扰动产生的高压金属流斜侵彻轨道的结果;随着速度的减小,冲击维持时间内峰值压力下降明显,如图 9所示,形成的金属流的动能也随之减小,对轨道的侵彻效果也变弱。
文献[6]中的实验也发现,在到达刨削临界速度之前,轨道上有磨损擦伤产生,如图 10所示,当速度超过0.8 km/s后比较明显,由磨损至产生刨削是一个逐渐演变的过程。文献[4-5]的轨道炮刨削实验,分别得出7075铝电枢和T2纯铜轨道的刨削临界速度为1.35和1.30 km/s,本文数值模拟的结果与之接近。由于刨削易发生在高速段,因此实际工程应用中,可通过提高高速段轨道的加工平整度或提高高速段轨道材料的强度来削弱及抑制刨削损伤的出现。
另外,本文中还计算分析了电枢头部倾角对刨削的影响,电枢速度及颗粒尺寸保持不变,改变电枢的头部倾角进行计算。结果显示,随着电枢倾角的增大,冲击形成的刨坑也变小变浅。这是因为,改变电枢头部倾角,相当于改变了电枢与微颗粒发生碰撞的角度,导致冲击碰撞所形成的高压金属流对轨道的侵彻角度发生变化,侵彻效果也自然受到影响。因此在工程实际应用中,设计合适的电枢头部倾角对削弱轨道的刨削损伤也很重要。
4. 结论
基于轨道炮结构特点及冲击热力学理论,采用物质点法构建了微颗粒诱发的刨削三维模型,对轨道刨削的形成过程进行了数值模拟,分析了刨削产生机理以及一些影响因素,得到以下主要结论:
(1) 刨削是在瞬态的局域高压扰动下产生的,电枢与轨道表面的局域高速冲击产生了瞬时的能量交换,使得受冲击的轨道材料变成高热高压的金属流斜角度侵彻下游轨道表面,塑性变形从冲击点开始迅速扩展,形成了泪滴状刨坑;
(2) 对于特定的枢轨材料和结构形式,刨削的产生存在速度阈值;超过该阈值,随着速度增加,刨削越严重;低于该阈值,轨道出现磨损擦伤,并随着冲击速度的进一步减小而减轻。另外,减小颗粒尺寸以及增大电枢头部倾角,均可降低刨削损伤;
(3) 基于本文中对轨道炮刨削的数值模拟结果,表明物质点法用于轨道刨削这类冲击动力学问题的计算具有突出的优点,一是便于处理高速滑移界面上的材料塑性流动大变形破坏问题,避免了有限元计算容易出现的网格畸变;二是在应变梯度很大的局域内不必划分很细密的网格也可以较方便地模拟出断裂失效等问题,有着较高的计算效率。
本文中,采用数值方法模拟了轨道刨削的形成过程并分析了颗粒大小、冲击速度、电枢倾角的影响,下一步的研究将更多地考虑枢轨材料与结构的影响。
期刊类型引用(4)
1. 杨小伟,王浩,王正宏. PBX-9407炸药冲击波比能流密度研究. 兵器装备工程学报. 2025(02): 130-135 . 百度学术
2. 杜宁,赵梓淇,熊玮,刘闯,张先锋. 壳体厚度对装药爆炸冲击波特性影响研究. 弹道学报. 2023(03): 72-77 . 百度学术
3. 高涵,冯晓军,尚宇,张坤. 混合炸药微结构设计与制备研究进展. 火炸药学报. 2023(09): 761-775 . 百度学术
4. 龚悦,汪旭光,何杰,颜事龙,程扬帆. 铝粉粒度对乳化炸药能量输出特性及热安定性的影响. 化工学报. 2017(04): 1721-1727 . 百度学术
其他类型引用(5)
-
计量
- 文章访问数: 2519
- HTML全文浏览量: 105
- PDF下载量: 358
- 被引次数: 9