一种长贮库房内偏二甲肼扩散浓度分布的快速估算方法与流程

本发明属于航天发射场安全性分析领域，是一种长贮库房内偏二甲肼扩散浓度分布的快速估算方法。

背景技术：

随着航天事业的发展，发射任务密集，发射场安全问题显得尤为重要。偏二甲肼为无色有毒易燃液体，是导弹、卫星、飞船等发射试验和运载火箭的主体燃料，根据任务需要，通常长期贮存一定量的推进剂。在偏二甲肼长贮库房内，管路连接处、阀门、法兰等位置存在发生泄漏的可能，泄漏后在地面蒸发为有毒气体进入空气中，对库房内工作人员的安全产生威胁。因此，有必要对泄漏事故进行快速评估，分析不同泄漏时间有毒推进剂在库房中的扩散范围以及人员危害区域，以便及时采取防护措施，保证人员的安全。

正确估算各种条件下危害性气体浓度随时空的分布和变化，关键在于扩散模型在各种条件下的应用。国内外关于不同条件下气体扩散模型的应用有比较多的研究，包括高架源气体泄漏、滨海地形大气扩散、城市污染气体扩散、火灾烟雾扩散等[任建国,鲁顺清.气体扩散数学模型在安全评价方面的应用[j].中国安全科学学报,2006,16(3):12-16.]。关于发射场推进剂泄漏，陈新华等提出了爆炸后残余推进剂蒸发形成的毒气在大气中扩散的解析解[陈新华,武江涛,佟连捷等.液体火箭爆炸后毒气扩散研究[j].推进技术,1999,20(5):6-10.]，但对于长贮库房内的扩散浓度计算方法未有相关报道。研究有毒气体扩散浓度最常用的方法是基于流体动力学的数值模拟与仿真，这种方法从建模到结果分析需要几个小时，甚至几天的时间，不适用于急需抢险救援的发射场推进剂长贮库房。

本发明的优势在于，可以根据长贮库房内偏二甲阱气体传感器报警数据，快速估算库房内偏二甲阱气体浓度分布，与基于流体动力学数值模拟获得长贮库房内偏二甲阱浓度分布的方法相比，本发明的一种长贮库房内偏二甲肼扩散浓度分布的快速估算方法，简单、易于实现、误差小、速度快。

技术实现要素：

本发明的一种长贮库房内偏二甲肼扩散浓度分布的快速估算方法，目的是根据传感器的报警数据，结合泄漏速率、平均风速、影响区域、地面反射系数和扩散参数等影响因素，快速估算偏二甲阱长贮库房内有毒气体的浓度分布，为制定快速安全的抢险堵漏方案提供有力支撑。

本发明提供了一种长贮库房内偏二甲肼扩散浓度分布的快速估算方法，包括：

步骤一、提出基本假设条件。泄漏物直接以气态或液相蒸发为气态的形式进入空气中，且在泄漏扩散过程中不存在相态转换；泄漏源为连续泄漏源，计算所得浓度为稳态浓度分布，不稳定流场造成的浓度波动不予考虑；扩散过程适用梯度输送理论，湍流引起的动量通量与风速梯度成正比，在平稳、均匀湍流下扩散浓度分布符合正态分布形式。

步骤二、确定偏二甲阱浓度分布基本计算公式。在偏二甲肼长贮库房中，主要是小流量的偏二甲阱渗漏，泄漏类型以点源和面源为主。在步骤一的基本假设条件下，根据气体扩散方程解析理论，分别建立点源和面源泄漏时，偏二甲阱浓度分布的基本计算公式。

步骤三、确定偏二甲阱浓度分布基本计算公式中的参数。偏二甲阱浓度分布基本计算公式中需要确定的参数包括泄漏源泄漏速率、平均风速、地面反射系数和扩散参数。

其特征在于：

(1)通过合理假设，使得所建偏二甲阱浓度分布基本计算公式简单、实用。

(2)充分考虑了储罐液面以上气相泄漏、储罐液面以下液相泄漏、阀门和管路破损发生小孔泄漏，以及阀门与管路或储罐的连接处等四种最有可能发生泄漏的情形对泄漏速率的影响。

(3)基于大量数值仿真数据，确定了平均风速计算模型。

(4)充分考虑长贮库房内不同区域对扩散参数的影响。

本发明的优势是：

可以根据偏二甲阱长贮库房内的偏二甲阱浓度传感器报警数据，快速估算库房内偏二甲阱的浓度分布，易于发射场工作人员操作使用，为制定抢险堵漏方案提供有力支撑。

附图说明

图1本发明的泄漏扩散浓度计算主要步骤示意图；

图2本发明的偏二甲肼长贮库房三维模型示意图；

图3本发明的储罐间速度分布曲线示意图；

图4本发明的多组风速拟合参数；

图5本发明的根据速度分布影响的区域划分示意图；

图6本发明的根据浓度分布影响的区域划分示意图；

图7本发明的不同下风向距离地面反射系数；

图8本发明的下风向不同距离地面反射系数及残差示意图；

图9本发明的不同稳定度的扩散参数经验公式；

图10本发明的采样点分布示意图；

图11本发明的泄漏位置下风向1.5m取样点浓度；

图12本发明的偏二甲肼浓度分布拟合曲线；

图13本发明的泄漏位置下风向0.5～7.5m的y方向扩散参数；

图14本发明的y方向扩散参数与下风向距离的拟合曲线及残差；

图15本发明的泄漏位置下风向0.5～7.5m的z方向扩散参数；

图16本发明的z方向扩散参数与下风向距离的拟合曲线及残差。

具体实施方式

结合附图对本发明的一种长贮库房内偏二甲肼扩散浓度分布的快速估算方法做进一步详细描述。

图1为本发明的泄漏扩散浓度计算主要步骤。

步骤一、提出基本假设条件。泄漏物直接以气态或液相蒸发为气态的形式进入空气中，且在泄漏扩散过程中不存在相态转换；泄漏源为连续泄漏源，计算所得浓度为稳态浓度分布，不稳定流场造成的浓度波动不予考虑；扩散过程适用梯度输送理论，湍流引起的动量通量与风速梯度成正比，在平稳、均匀湍流下扩散浓度分布符合正态分布形式。

步骤二、确定偏二甲阱浓度分布基本计算公式。在偏二甲肼长贮库房中，主要是小流量的偏二甲阱渗漏，泄漏类型以点源和面源为主。在步骤一的基本假设条件下，根据气体扩散方程解析理论，分别建立点源和面源泄漏时，偏二甲阱浓度分布的基本计算公式。

以泄漏源为原点，下风向为x轴建立三维坐标系，如图2所示，x，y，z为所求位置坐标；点源泄漏指泄漏孔径较小，泄漏孔径对泄漏扩散的影响可以忽略不计，泄漏扩散浓度基本计算公式为

式中，t为泄漏时间；q为泄漏量，与泄漏时间相关；u为泄漏位置风速；σy和σz分别为y和z方向的扩散参数，λ为地面反射系数。

面源泄漏空间浓度分布计算公式为

式中，t为泄漏时间；q为泄漏量；u为泄漏位置风速；σy和σz分别为y和z方向的扩散参数，λ为地面反射系数。

步骤三、确定偏二甲阱浓度分布基本计算公式中的参数。偏二甲阱浓度分布基本计算公式中需要确定的参数包括泄漏源泄漏速率、平均风速、地面反射系数和扩散参数。

(1)泄漏速率计算

泄漏速率计算采用成熟经验公式，主要包括如下几种状态：

1)储罐液面以上气相泄漏，由于孔洞自由泄放的气体泄漏量计算公为

式中，p为储罐内绝对压力，pa；r为理想气体常数，8.314j/(mol·k)；m为气体摩尔质量，kg/mol；ρ为密度，kg/m³；t为温度，k；γ为绝热指数。c0为泄漏系数，一般取值范围为0.6～1.0。

2)储罐液面以下液相泄漏时，若不发生闪蒸的泄漏量的计算公式为

式中，cd为排放系数，取决于孔的形状和流动状态，液体一般为0.6～0.64；ar为泄漏孔的有效面积，m²；ρ为有害液体的密度，kg/m³；g为重力加速度，m/s²；h为泄漏位置距液面的高度，m；p1为容器内部压力，pa；pa为外界大气压强，pa。

3)当阀门、管路破损发生小孔泄漏时，泄漏强度的计算方法为

式中，q为泄漏量，kg/s；a为泄漏口面积，m²；c0为扩散系数，ρ为液体密度，kg/m³；p1和p2分别表示泄漏位置上流压强和下流压强，pa。

机械或外力方面的原因也可能造成阀门故障，如泄压元件安全阀失效导致故障排气的发生。此时气体的泄漏量计算公式为

式中，c0为流量系数，与阀门排气口结构有关；x为气体特征系数；av为阀门排气口面积，mm²；m为泄漏物质的摩尔质量，kg/kmol；z为气体的压缩系数；t为气体温度，k。

4)密封件或连接件的失效导致阀门与管路或储罐的连接处以及管道连接处发生泄漏，泄漏量计算公式为

式中，q为泄漏量，kg/s；h为间隙的大小，cm；δp为设备内外的压力差，mpa；ρ为密度，kg/m³；η为介质的粘度，pa·s；r1为内径，cm；r2为外径，cm。

(2)平均风速计算

偏二甲肼长贮库房在对称通风条件下风速分布规律明显，可通过数值模拟结果得到风速的分布规律和计算方法，以公式的形式纳入计算。

按照进风口的中心线风速衰减规律，风速应满足公式

式中，u为风速，u0为初始风速，s为与进风口之间距离，a1、b1为常系数，a为声速。

取10组储罐通风口中心线0～16.4m速度分布，排除排风口附近2m对风速影响，按照风速衰减规律进行拟合如图3所示。

经过多组数据的拟合，得到拟合参数值如图4所示，取均值因此，在该通风条件下，平均风速的计算可以通过该计算公式近似计算。

式中，u0＝4.306m/s，a＝340m/s，s为与进风口的距离。

(3)影响区域划分

按照壁面与储罐对推进剂扩散规律影响研究的结果，不同区域的浓度分布规律是不同的，其气体扩散浓度计算公式的参数也不同，故需要对不同区域进行划分，根据数值模拟结果分别对浓度分布的计算方法进行分析。

壁面及储罐主要从空间上对推进剂扩散过程起到阻碍作用，影响局部区域流速的分布及大小，进而对扩散浓度的分布规律产生影响。图5为根据速度分布的影响区域划分，横向主要按照进风口和排风口位置划分10个部分，标号为ⅰ～ⅹ，其中，ⅰ～ⅸ区域的速度分布基本一致，ⅹ区域的速度分布受到地下空间的影响，速度相对较小，需单独考虑分析其速度分布。

图6为根据浓度分布的影响区域划分，不同泄漏位置的浓度分布规律是相似的，以第二与第三储罐之间泄漏为例，主要划分为3个区域，在图6中标记为①～③。

①区域主要以自由扩散为主，由于泄漏位置距离进风口较近，y方向浓度分布范围较小，x方向主要向下风向扩散，浓度分布规律与理想气体扩散方程所得的浓度分布形式较为接近，可利用理想气体扩散方程的修正形式较为准确的计算该区域的浓度分布。②区域则受到储罐及支架的阻碍作用，扩散浓度分布呈收缩状，需对扩散公式进行较大的收缩性修正。③区域受到储罐扰流，壁面阻挡作用及排风口的排出作用，浓度分布不稳定，较低浓度值分布范围较大，中心浓度则在排风口附近收缩，存在随时间左右摆动的现象，有一定不稳定性。

(4)地面反射系数计算

沿下风向每隔0.5米取50取样点，最大下风向距离为9m，分别求出地面反射系数，如图7所示。

对地面反射系数与下风向距离进行拟合，得到的拟合曲线如图8(a)所示，拟合曲线方程为

式中，y0＝0.01535，a2＝2.02934，w1＝2.04311，xc＝1.89164，拟合残差如图8(b)所示。

(5)扩散参数的确定与修正

长贮库房内气体流动基本保持不变，按照湍流稳定度的分类方法，长贮库房内大部分区域平均风速在1～2m左右，室内无太阳辐照，湍流稳定度基本在c类和d类之间。开放空间的扩散参数经验公式如图9所示，c类和d类的表达形式相近，区别在于经验公式中的系数。开放空间的扩散参数经验公式不能直接适用于长贮库房内的推进剂泄漏扩散计算，但由于湍流扩散的本质是相同的，扩散参数的表达形式应是相同的，可根据长贮库房内浓度分布规律对开放空间扩散经验公式中的参数进行一定修正，以适应长贮库房内扩散分布的实际情况。

y方向和z方向扩散参数表达式为

σy＝a1x/(1+b1x)^0.5，σz＝a2x/(1+b2x)^0.5(11)

为了得到长贮库房内偏二甲阱扩散参数的表达式，分别对沿下风向不同位置的地面浓度分布及垂直浓度分布进行取样和高斯分布拟合，如图10所示，y方向和z方向的扩散参数求解方法主要取下风向0.5m～7.5m的y方向地面浓度分布数据进行曲线拟合。

以泄漏下风向1.5m的y方向地面浓度分布数据为例，取地面距离泄漏位置1.5m处，宽度为4m的范围内，采集取样点的浓度数据，如图11所示，共设置21个取样点，相邻两个取样点之间的距离为0.2m。按照湍流扩散和扩散统计理论，浓度在y方向呈正态分布，其特征值σ代表了气体在y方向的扩散能力，与公式中扩散参数σy的物理意义是相同的，其值相等。

如图12，对取样点浓度进行高斯分布拟合，两侧浓度偏高是由于两侧位置位于储罐底部，受到储罐的影响，导致扩散空间有限，故浓度偏高。经过多次迭代，拟合曲线所得表达式为

式中，w2＝2σ＝0.66665，an＝0.00144，σ＝0.33332。

将x＝1.5m，z＝0m，t＝300s，代入公式(1)计算得

求得与an相应的系数为0.0014097，这与拟合结果基本一致，间接验证了公式的正确性。根据曲线拟合结果，得到泄漏位置下风向1.5m处y方向扩散参数σ为0.33332。

采用同样的处理方法，对下风向0.5m～7.5m范围内y方向地面浓度分布数据进行曲线拟合，结果如图13所示。

以d类稳定度的扩散参数经验公式形式对数值模拟所得长贮库房内推进剂泄漏扩散的扩散参数进行拟合，拟合方程为

σy＝m1δx/(1+n1δx)^0.5(14)

式中，m1和n1为待拟合参数，σy为y方向扩散参数，x为泄漏位置的下风向距离。经过拟合，得到参数m1为165.37，n1为376383.16。表征拟合度的残差r²为0.99079，超过0.99，如图14所示，拟合程度基本满足要求。因此，长贮库房内y方向偏二甲阱扩散参数与下风向距离的关系式为

σy＝165.37δx/(1+376383.16δx)^0.5(15)

利用同样的方法对z方向扩散参数进行拟合求解，通过拟合得到的z方向扩散参数值如图15所示。

对z方向扩散参数和下风向距离x进行拟合，拟合函数为

σz＝m2δx/(1+n2δx)^0.5(16)

式中，m2和n2为待拟合参数，σz为z方向扩散参数，x为泄漏位置的下风向距离。经过拟合，得到参数m2为0.2912，n2为2.2597。表征拟合度的残差r²为0.98687，接近0.99，如图16所示，拟合程度基本满足要求。

因此，长贮库房内z方向扩散参数与下风向距离的关系式为

σz＝0.2912δx/(1+2.2597δx)^0.5(17)

以上所得扩散参数是根据推进剂长贮库房内泄漏扩散的数值模拟结果拟合所得，某一位置的扩散参数与该位置距泄漏源位置的下风向距离有关。由于不同位置储罐对其影响不同，在储罐间误差有限，该结果可以应用到储罐间不同位置的泄漏扩散浓度估算中。