CN104834825A - 一种卧式柱形爆炸源碎片击中概率评估方法 - Google Patents

Abstract

本发明公开了一种卧式柱形爆炸源碎片击中概率评估方法，其是基于Monte-Carlo模拟的方法，考虑卧式柱形爆炸源尺寸和支架高度，从统计意义上计算爆炸碎片的击中其它目标出储罐的概率。本发明提出的击中概率评估方法考虑卧式柱形爆炸源尺寸和支架高度，具有较高的计算精度，填补了柱形容器这一复杂形状爆炸源击中概率计算方法的空白。

Description

一种卧式柱形爆炸源碎片击中概率评估方法

技术领域

本发明涉及石油化工技术领域，具体涉及一种导致多米诺效应的评估方法，特别是由爆炸碎片击中其他目标概率的计算方法。

背景技术

液化石油气(LPG)、液化天然气(LNG)等带压储罐在火灾、撞击等外部能量作用下，可能会导致沸腾液体扩展蒸气爆炸(BLEVE)，BLEVE的发生常伴随热辐射、爆炸冲击波和碎片，而热辐射、爆炸冲击波和碎片则是造成多米诺效应最主要的原因，尤其是爆炸产生的具有一定初始速度的碎片，在飞行过程中容易可能击中人和重要设备，进而导致更大的灾难性后果。

如何评估碎片击中其他目标的概率，发展出了二维方程和三维方程模型的方法，使得计算获得的击中概率越来越精确。

“A procedure for analyzing the flight of missiles from explosions ofcylindrical vessels”基于二维方程模型，对碎片抛射的速度、质量等随机性进行了分析，探讨了由圆柱形爆炸源产生的碎片不同抛射距离范围内的落地概率。“Integrated probabilistic framework for domino effect and risk analysis”针对爆炸碎片击中目标设备概率的各种计算方法，提出了击中概率三维方程模型，即当爆炸碎片的抛射轨迹和目标物体方程相交时，二者发生碰撞。“爆炸球罐尺寸对抛射碎片击中相邻罐体概率的影响”通过将球罐的尺寸考虑在内，建立了新三维的击中概率计算方法。此外，通过以碎片、柱形和概率为关键词在国家知识产权局网站对国内外专利的进行索，并无发现类似或相关的专利技术。

然而，国外已有的研究都将爆炸物体当作质点考虑，而国内已有的研究也都未给出考虑卧式柱形爆炸源尺寸的击中概率的计算方法，并且也都忽略了储罐支撑架高度的影响。

本发明就是为解决上述存在的技术问题而提出的。

发明内容

为解决现有技术中存在的问题，本发明提出的击中概率评估方法考虑卧式柱形爆炸源尺寸和支架高度，填补了柱形容器这一复杂形状爆炸源击中概率计算方法的空白，并将储罐支撑架的高度纳入到计算方法中。

本发明使得击中概率计算结果更加精确、更加符合实际情况，更可以为评估装置的危险性、合理布局及采取相应的保护措施提供准确的依据。

本发明提供一种卧式柱形爆炸源碎片击中概率评估方法，其是基于Monte-Carlo模拟的方法，考虑卧式柱形爆炸源尺寸和支架高度，从统计意义上计算爆炸碎片的击中概率。

其中，所述评估方法包括：

S1步：输入爆炸源的基本参数，包括卧式柱形爆炸罐体积V，半径R、长度L、总质量M、材料密度ρ以及空气密度ρ₀和大气压力P₀；

S2步：爆炸模拟次数n进行递增，n初始值为1；

S3步：判断爆炸模拟次数n是否大于设定最大模拟次数N，若是则执行步骤S9；否则执行步骤S4；

S4步：通过随机数发生器，产生随机种子，计算爆炸产生的碎片初始参数，包括碎片的水平抛射角θ、垂直抛射角碎片初速度v₀、阻力系数K_D和升力系数K_L；

S5步：将S4步获得的已有的随机参数作为初始值，结合碎片轨迹三维二阶导数方程计算获得碎片的轨迹方程；

S6步：假设目标的中心坐标为(x₁，0，0)，那么对于目标是球罐的体积方程是(x-x₁)²+y²+z²＝r²；

S7：联立S5步和S6步所获得的方程，采用二分法，求解方程组是否有解，若方程有解，执行S8步，反之则执行S2步；

S8：增加一次到爆炸碎片击中目标的次数N1上；

S9：爆炸碎片击中目标的概率P＝N1/N，流程结束。

其中，所述S4步进一步具体包括：

S4-1步：计算产生的碎片水平抛射角θ，其服从20％-[30，150°]；30％-[150，210°]；20％-[210，330°]；30％-[330，30°]的均匀分布；

S4-2步：计算产生的碎片垂直抛射角其是[-90，90°]上的

S4-3步：计算产生的碎片产生的数量n，分布函数是其中λ₀＝0.8145，λ₁＝0.2252，λ₃＝0.0321；

S4-4步：计算爆炸时的失效压力P₁，值为操作压力的[0.9，1.1]倍上均匀分布；

S4-5步：计算爆炸转化到碎片上的能量系数c，函数分布是[0.2，0.5]上的右三角分布；

S4-6步：爆炸产生各个碎片质量m，碎片质量服从分布，其中a＝0.41213，

S4-7步：计算阻力系数C_L和升力系数C_D，其中，若碎片为端盖和长端盖则C_D服从[0.8，1.1]上的均匀分布，若碎片是平板形则C_D应服从[1.1，1.8]上的均匀分布；而对于碎片是端盖和长端盖形状C_L则服从[0.351，0.468]均匀分布，对于是平板时，则C_L为0；

S4-8步：通过S4-4步和S4-5步获得的失效压力P₁和能量比例系数c，爆炸产生的能量 $E=\[1-{\left(\frac{p_0}{p_1}\right)}^{(\mathrm{\γ}-1)/\mathrm{\γ}}+(\mathrm{\γ}-1)\frac{p_0}{p_1}\]\frac{p_1}{\mathrm{\γ}-1}V,$

其中：P₀是爆炸环境的大气压力；γ＝1.33，是比热系数，而后计算爆炸碎片从爆炸中获得的能量E_k＝c×E；

S4-9步：由碎片的能量E_k、碎片数量n，可以获得碎片所获得的动能；

S4-10步：依据碎片的动能和质量m，可以获得碎片的初速度

S4-11步：由S4-7步获得的C_L、C_D值，可以计算二阶导数方程中的阻力参数 $k_D=\frac{{\mathrm{\ρC}}_D{A}_D}{2m},$ 升力参数 $k_L=\frac{{\mathrm{\ρC}}_L{A}_L}{2m}.$

其中，所述S5步进一步包括：

S5-1步：获取S4-10步计算获得的初始参数，碎片初速度v₀；

S5-2步：获取S4-11步计算获得的初始参数二阶导数方程中的阻力系数k_D和升力系数k_L；

S5-3步：输入作为储罐底座的支撑架高度h；

S5-4步：获取爆炸源尺寸参数，包括半径R和柱形罐长度L；

S5-5步：获取S4-2步计算获得的水平抛射角θ；

S5-6步：获取S4-3步计算获得的随机垂直抛射角

S5-7步：根据S5-3步至S5-6步所获得的参数值，根据爆炸源在坐标系中的具体位置，计算其中心点的位置(x₀，y₀，z₀)；

S5-8步：根据S5-7步和S5-1步、S5-2步获得的爆炸碎片初始位置和碎片初速度、阻力系数等，参照二阶导数方程 $\left\{\begin{array}{c}\frac{d^{2}x}{{\mathrm{dt}}^2}+\left(k_D-{\left(-1\right)}^q{k}_L\right){\left(\frac{dx}{dt}\right)}^2=0\\ \frac{d^{2}z}{{\mathrm{dt}}^2}+\left(k_D-{\left(-1\right)}^q{k}_L\right){\left(\frac{dz}{dt}\right)}^2=0\\ \frac{d^{2}y}{{\mathrm{dt}}^2}+\left({\left(-1\right)}^q{k}_D-k_L\right){\left(\frac{dy}{dt}\right)}^2+g=0\end{array}\right.,$ 进行积分求解，获取爆炸碎片抛射的轨迹方程。

其中，所述S7步进一步具体包括：

S7-1步：以爆炸发生时刻，作为零时刻起点，设定一个极小的时间步作为爆炸碎片正在飞行的时刻；

S7-2步：将碎片飞行时刻作为输入初始值，代入到S5步所获得的碎片轨迹方程中，获得碎片该时刻位置(x，y，z)；

S7-3步：计算上一时刻与该时刻的时间差值是否满足精度要求，若是则执行S7-4步，反之，则执行S7-5步；

S7-4步：流程结束，方程无解；

S7-5步：根据碎片所处位置(x，y，z)，判断其是否在目标储罐的内部，若是，则执行S7-6步，反之，执行S7-7步；

S7-6步：方程有解，流程结束；

S7-7步：利用二分法，对碎片飞行时刻进行更新，而后执行S7-2步。

本发明有益的技术效果在于：

本发明提出的卧式柱形爆炸源碎片集中概率评估方法使得击中概率计算结果更加精确，尤其在评估爆炸源附近装置的危险性时，可给予更为准确的指导，进而使得装置布局合理、保护措施具有准确的针对性。

附图说明

图1为本发明集中概率评估方法的计算流程图；

图2为基于Monte Carlo碎片随机参数获取；

图3为碎片轨迹方程求解流程；

图4为方程组联立求解。

具体实施方式

本发明提供一种卧式柱形爆炸源碎片击中概率评估方法，是基于Monte-Carlo模拟的方法，考虑卧式柱形爆炸源尺寸和支架高度，从统计意义上计算爆炸碎片的击中概率，以下结合附图和实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施(图1所示)。

S1步：输入爆炸源的基本参数，包括卧式柱形爆炸罐体积V＝130m³，半径R、长度L、总质量M、材料密度ρ以及空气密度ρ₀和大气压力P₀＝101325Pa；

S2步：爆炸模拟次数n进行递增，起初n＝1；

S3步：判断爆炸模拟次数n是否大于设定最大模拟次数N＝200000，若是则执行步骤S9步；否则执行步骤S4步；

S4步：通过随机数发生器，产生随机种子，计算爆炸产生的碎片初始参数，包括碎片的水平抛射角θ、垂直抛射角碎片初速度v₀、阻力系数K_D和升力系数K_L，具体步骤执行S4-1步至S4-11步；

S4-1步：计算产生的碎片水平抛射角θ，其服从20％-[30，150°]；30％-[150，210°]；20％-[210，330°]；30％-[330，30°]的均匀分布；

S4-2步：计算产生的碎片垂直抛射角其是[-90，90°]上的

S4-3步：计算产生的碎片产生的数量n，分布函数是其中λ₀＝0.8145，λ₁＝0.2252，λ₃＝0.0321；

S4-4步：计算爆炸时的失效压力P₁，值为操作压力的[0.9，1.1]倍上均匀分布；

S4-5步：计算爆炸转化到碎片上的能量系数c，函数分布是[0.2，0.5]上的右三角分布；

S4-6步：爆炸产生各个碎片质量m，碎片质量服从 $p_k\left(k\right)=1_{\[0,1\]}\left(k\right)\frac{1}{B\left(a,b\right)}{k}^{a-1}{\left(1-k\right)}^{b-1}$ 分布，其中a＝0.41213，

S4-7步：计算阻力系数C_L和升力系数C_D，其中，若碎片为端盖和长端盖则C_D服从[0.8，1.1]上的均匀分布，若碎片是平板形则C_D应服从[1.1，1.8]上的均匀分布；而对于碎片是端盖和长端盖形状C_L服从[0.351，0.468]均匀分布，对于是平板时，则C_L为0；

S4-8步：通过S4-4步和S4-5步获得的失效压力P₁和能量比例系数c，爆炸产生的能量 $E=\[1-{\left(\frac{p_0}{p_1}\right)}^{(\mathrm{\γ}-1)/\mathrm{\γ}}+(\mathrm{\γ}-1)\frac{p_0}{p_1}\]\frac{p_1}{\mathrm{\γ}-1}V,$

其中：P₀是爆炸环境的大气压力；γ＝1.33，是比热系数，

而后计算爆炸碎片从爆炸中获得的能量E_k＝c×E；

S4-9步：由碎片的能量E_k、碎片数量n，可以获得碎片所获得的动能；

S4-10步：依据碎片的动能和质量m，可以获得碎片的初速度

S4-11步：由S4-7步获得的C_L、C_D值，可以计算二阶导数方程中的阻力参数 $k_D=\frac{{\mathrm{\ρC}}_D{A}_D}{2m},$ 升力参数 $k_L=\frac{{\mathrm{\ρC}}_L{A}_L}{2m};$

S5步：将S4步获得的已有的随机参数作为初始值，结合碎片轨迹三维二阶导数方程计算获得碎片的轨迹方程，具体步骤执行S5-1步至S5-8步；

S5-1步：获取S4-10步计算获得的初始参数，碎片初速度v₀；

S5-2步：获取S4-11步计算获得的初始参数二阶导数方程中的阻力系数k_D和升力系数k_L；

S5-3步：输入作为储罐底座的支撑架高度h；

S5-4步：获取爆炸源尺寸参数，包括半径R和柱形罐长度L；

S5-5步：获取S4-2步计算获得的水平抛射角θ；

S5-6步：获取S4-3步计算获得的随机垂直抛射角

S5-7步：根据S5-3步至S5-6步所获得的参数值，根据爆炸源在坐标系中的具体位置，这里假设爆炸源x、y方向坐标都为0，然后计算其碎片初始的位置(x₀，y₀，z₀)，其中垂直于轴向坐标轴向坐标z₀＝l；纵向坐标中间参数其中，l在[-L/2，L/2]均匀分布；

S5-8步：根据S5-7步和S5-1步、S5-2步获得的爆炸碎片初始位置和碎片初速度、阻力系数等，参照二阶导数方程 $\left\{\begin{array}{c}\frac{d^{2}x}{{\mathrm{dt}}^2}+\left(k_D-{\left(-1\right)}^q{k}_L\right){\left(\frac{dx}{dt}\right)}^2=0\\ \frac{d^{2}z}{{\mathrm{dt}}^2}+\left(k_D-{\left(-1\right)}^q{k}_L\right){\left(\frac{dz}{dt}\right)}^2=0\\ \frac{d^{2}y}{{\mathrm{dt}}^2}+\left({\left(-1\right)}^q{k}_D-k_L\right){\left(\frac{dy}{dt}\right)}^2+g=0\end{array}\right.,$ 进行积分求解，获取爆炸碎片抛射的轨迹方程(如表1)；

S6步：假设目标为球罐，体积为1600m³，则目标中心坐标为(x₁，0，0)，目标方程是(x-x₁)²+y²+z²＝r²；

S7步：联立S5步和S6步所获得的方程，采用二分法，求解方程组是否有解，具体执行步骤S7-1步至S7-7步。若方程有解，执行S8步，反之则执行S2步；

S7-1步：以爆炸发生时刻，作为零时刻起点，设定一个极小的时间步0.0005s作为爆炸碎片初始正在飞行的时刻；

S7-2步：将碎片飞行时刻作为输入初始值，代入到S5步所获得的碎片轨迹方程中，获得碎片该时刻位置(x，y，z)；

S7-3步：计算上一时刻与该时刻的时间差值Δt是否满足精度0.001s要求，若是则执行S7-4步，反之，则执行S7-5步；

S7-4步：流程结束，方程无解；

S7-5步：根据碎片所处位置(x，y，z)，判断其是否在目标储罐的内部，若是，则执行S7-6步，反之，执行S7-7步；

S7-6步：方程有解，流程结束；

S7-7步：利用二分法重新计算碎片飞行时刻，判断该时刻碎片位置x＜x₁若是t_low＝t，t＝(t_high+t_low)/2，反之t_high＝t，t＝(t_high+t_low)/2，其中，t_high为碎片飞行时刻上限范围值，t_low为碎片飞行时刻下限范围值，而后执行S7-2步；

S8步：增加一次到爆炸碎片击中目标的次数N1上；

S9步：爆炸碎片击中目标的概率P＝N1/N。

表1碎片轨迹方程

其中，km＝k_D-k_L；kp＝k_D+k_L

u_x/u_y/u_z是x/y/z方向碎片的初速度；

x_p/y_p/z_p是达到最高点的坐标；

tp碎片达到最高点的时间。

所有上述的首要实施这一知识产权，并没有设定限制其他形式的实施这种新产品和/或新方法。本领域技术人员将利用这一重要信息，上述内容修改，以实现类似的执行情况。但是，所有修改或改造基于本发明新产品属于保留的权利。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例。但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。

Claims (5)

1.一种卧式柱形爆炸源碎片击中概率评估方法，其特征在于：是基于Monte-Carlo模拟的方法，考虑卧式柱形爆炸源尺寸和支架高度，从统计意义上计算爆炸碎片的击中概率。

2.如权利要求1所述的概率评估方法，其特征在于，包括：

S1步：输入爆炸源的基本参数，包括卧式柱形爆炸罐体积V，半径R、长度L、总质量M、材料密度ρ以及空气密度ρ₀和大气压力P₀；

S2步：爆炸模拟次数n进行递增，n初始值为1；

S3步：判断爆炸模拟次数n是否大于设定最大模拟次数N，若是则执行步骤S9；否则执行步骤S4；

S4步：通过随机数发生器，产生随机种子，计算爆炸产生的碎片初始参数，包括碎片的水平抛射角θ、垂直抛射角碎片初速度v₀、阻力系数K_D和升力系数K_L；

S5步：将S4步获得的已有的随机参数作为初始值，结合碎片轨迹三维二阶导数方程计算获得碎片的轨迹方程；

S6步：假设目标的中心坐标为(x₁，0，0)，那么对于目标是球罐的体积方程是(x-x₁)²+y²+z²＝r²；

S7：联立S5步和S6步所获得的方程，采用二分法，求解方程组是否有解，若方程有解，执行S8步，反之则执行S2步；

S8：增加一次到爆炸碎片击中目标的次数N1上；

S9：爆炸碎片击中目标的概率P＝N1/N，流程结束。

3.如权利要求1或2所述的概率评估方法，其特征在于：所述S4步进一步具体包括，

S4-1步：计算产生的碎片水平抛射角θ，其服从20％-[30，150°]；30％-[150，210°]；20％-[210，330°]；30％-[330，30°]的均匀分布；

S4-2步：计算产生的碎片垂直抛射角其是[-90，90°]上的

S4-3步：计算产生的碎片产生的数量n，分布函数是其中λ₀＝0.8145，λ₁＝0.2252，λ₃＝0.0321；

S4-4步：计算爆炸时的失效压力P₁，值为操作压力的[0.9，1.1]倍上均匀分布；

S4-5步：计算爆炸转化到碎片上的能量系数c，函数分布是[0.2，0.5]上的右三角分布；

S4-6步：爆炸产生各个碎片质量m，碎片质量服从分布，其中a＝0.41213，

S4-7步：计算阻力系数C_L和升力系数C_D，其中，若碎片为端盖和长端盖则C_D服从[0.8，1.1]上的均匀分布，若碎片是平板形则C_D应服从[1.1，1.8]上的均匀分布；而对于碎片是端盖和长端盖形状C_L则服从[0.351，0.468]均匀分布，对于是平板时，则C_L为0；

S4-8步：通过S4-4步和S4-5步获得的失效压力P₁和能量比例系数c，爆炸产生的能量 $E=[1-{\left(\frac{p_0}{p_1}\right)}^{(\mathrm{\γ}-1)/\mathrm{\γ}}+(\mathrm{\γ}-1)\frac{p_0}{p_1}]\frac{p_1}{\mathrm{\γ}-1}V,$

其中：P₀是爆炸环境的大气压力；γ＝1.33，是比热系数，而后计算爆炸碎片从爆炸中获得的能量E_k＝c×E；

S4-9步：由碎片的能量E_k、碎片数量n，可以获得碎片所获得的动能；

S4-10步：依据碎片的动能和质量m，可以获得碎片的初速度

S4-11步：由S4-7步获得的C_L、C_D值，可以计算二阶导数方程中的阻力参数 $k_D=\frac{\mathrm{\ρ}{C}_D{A}_D}{2m},$ 升力参数 $k_L=\frac{\mathrm{\ρ}{C}_L{A}_L}{2m}.$

4.如权利要求1至3所述的概率评估方法，其特征在于：所述S5步进一步包括，

S5-1步：获取S4-10步计算获得的初始参数，碎片初速度v₀；

S5-2步：获取S4-11步计算获得的初始参数二阶导数方程中的阻力系数k_D和升力系数k_L；

S5-3步：输入作为储罐底座的支撑架高度h；

S5-4步：获取爆炸源尺寸参数，包括半径R和柱形罐长度L；

S5-5步：获取S4-2步计算获得的水平抛射角θ；

S5-6步：获取S4-3步计算获得的随机垂直抛射角

S5-7步：根据S5-3步至S5-6步所获得的参数值，根据爆炸源在坐标系中的具体位置，计算其中心点的位置(x₀，y₀，z₀)；

S5-8步：根据S5-7步和S5-1步、S5-2步获得的爆炸碎片初始位置和碎片初速度、阻力系数等，参照二阶导数方程 $\left\{\begin{array}{c}\frac{d^{2}x}{{\mathrm{dt}}^2}+(k_D-{(-1)}^q{k}_L){\left(\frac{\mathrm{dx}}{\mathrm{dt}}\right)}^2=0\\ \frac{d^{2}z}{{\mathrm{dt}}^2}+(k_D-{(-1)}^q{k}_L){\left(\frac{\mathrm{dz}}{\mathrm{dt}}\right)}^2=0\\ \frac{d^{2}y}{{\mathrm{dt}}^2}+({(-1)}^q{k}_D-k_L){\left(\frac{\mathrm{dy}}{\mathrm{dt}}\right)}^2+g=0\end{array}\right.,$ 进行积分求解，获取爆炸碎片抛射的轨迹方程。

5.如权利要求1至4所述的概率评估方法，其特征在于：所述S7步进一步具体包括，

S7-1步：以爆炸发生时刻，作为零时刻起点，设定一个极小的时间步作为爆炸碎片正在飞行的时刻；

S7-2步：将碎片飞行时刻作为输入初始值，代入到S5步所获得的碎片轨迹方程中，获得碎片该时刻位置(x，y，z)；

S7-3步：计算上一时刻与该时刻的时间差值是否满足精度要求，若是则执行S7-4步，反之，则执行S7-5步；

S7-4步：流程结束，方程无解；

S7-5步：根据碎片所处位置(x，y，z)，判断其是否在目标储罐的内部，若是，则执行S7-6步，反之，执行S7-7步；

S7-6步：方程有解，流程结束；

S7-7步：利用二分法，对碎片飞行时刻进行更新，而后执行S7-2步。