CN104182588B - 一种重气持续泄漏扩散的动态模拟方法 - Google Patents

Abstract

本发明提供一种重气持续泄漏扩散的动态模拟方法，其是基于平板模型的，包括：输入已知参数，并根据已知参数初始化相关参数；计算ERPGs浓度稳态下的最大影响距离；在本地坐标系下生成左边界，左边界包括上风向左边界、下风向左边界和末端左边界；并校验纠正左边界的边界点；根据计算得到的左边界的边界点生成右边界；将左边界和右边界的边界点的坐标转换为经纬度坐标；并根据该经纬度坐标绘制出重气持续泄漏扩散的动态模拟GIS效果图。通过本发明，不仅仅可以模拟出稳定状态下的重气持续泄漏扩散态势，还可以动态的模拟重气持续泄漏扩散随时间变化的过程。

Description

一种重气持续泄漏扩散的动态模拟方法

技术领域

本发明涉及一种突发性事故的应急处理技术，特别是涉及一种重气持续泄漏扩散的动态模拟方法。

背景技术

危化品泄漏事故一般突发性强，影响范围广，危害后果严重；尤其是重气泄露事故，如氯气等，由于密度比空气大，泄漏后重气先重力下沉而后沿地面扩散，再加上泄漏物质本身的毒性，将对泄漏影响区域的人员、环境安全造成极大的威胁，甚至有可能造成大量的人员伤亡。因此，泄漏事故的控制越来越受社会的重视，国内外研究人员针对重气持续泄漏扩散也先后展开了一系列的研究，如现场试验，风洞试验，数学模拟等；相比于前两者方法，数值模拟的方法凭着其成本低，操作方便等优点以及计算机仿真技术的快速发展，近年来逐渐成为重气持续泄漏扩散研究的主流方法。

数学模拟法是指利用数学建模的方式来描述气体扩散过程，并利用计算机仿真技术实现对气体扩散态势的模拟；该方法大体分为两个步骤：1.构建重气扩散模型；2根据模型通过计算机仿真模拟气体扩散的过程。其中，扩散模型的构建起着至关重要的作用，它的准确性和科学性将直接决定最终模拟出来的结果，所以大量的国内外相关机构与学者围绕着构建更可靠更精确的扩散模型展开了大量的研究，目前，已经研究出了大量的重气扩散模型，如BM模型(经验模型)、箱模型、平板模型、浅层模型、拉格朗日粒子模型、CFD(Computational Fluid Dynamics，计算流体动力学)模型等。这些模型中，较复杂的模型，如CFD、浅层模型等，具有较好的准确性，然而由于其计算复杂度过大，计算时间过长，在实际的工业应用中，并不适用；反而是箱模型、平板模型等这样一些能够更好的兼顾计算复杂度和准确度的模型，更适合实际的工业应用。尤其是在对时间要求比较高的突发事故的应急处置中，平板模型的应用比较广泛。

平板模型是用于描述重气持续泄漏场景的重气持续泄漏扩散模型，其在模型的准确性和计算复杂度上具有很好的折中，很适合实际应用；然而由于该模型自身存在的一些数学上的限制性：

1、该模型属于是下风向方向上的一维平均分布模型，无法区分出垂风向上不同浓度的分布；

2、该模型是半封闭的模型，没有准确标识区域间的边界；

3、该模型是对重气持续泄漏扩散稳态的描述，不具备随时间变化的扩散过程的信息。

因此，基于平板模型使用计算机仿真重气持续泄漏扩散具有一定的难度，如采用一般的方法，简单、机械的对模型进行求解，仿真效果往往不够准确和直观，导致对突发事故应急处置的指挥调度、救援等工作的指导意义不大。

发明内容

鉴于以上所述现有技术的缺点，本发明的目的在于提供一种重气持续泄漏扩散的动态模拟方法，用于解决现有技术中通过平板模型无法动态模拟出随时间变化的重气持续泄漏扩散过程以及仿真结构不够准确和直观的问题。

为实现上述目的及其他相关目的，本发明提供一种重气持续泄漏扩散的动态模拟方法，所述方法是基于平板模型的，其特征在于，所述动态模拟方法包括：

步骤S10，输入已知参数，并根据所述已知参数初始化相关参数；

步骤S20，计算ERPGs浓度稳态下的最大影响距离；

步骤S30，在本地坐标系下生成左边界；

步骤S40，校验纠正所述左边界；

步骤S50，生成右边界：右边界的边界点的x坐标的值与所述左边界的边界点的x坐标的值相同，y坐标的值相反；

步骤S60，转换坐标系：将所述左边界和所述右边界的边界点的坐标转换为经纬度坐标；

步骤S70，输出所述步骤S60转换得到的经纬度坐标，并据此绘制GIS效果图。

可选地，所述已知参数包括风速、风向、气压、气温、泄漏源的经纬度、泄漏源的高度、重气云团初始密度、重气云团初始半宽b₀、泄漏后的时间、重气的分子量和ERPGs浓度。

可选地，所述相关参数包括：用弧度表示的风向；气压与标准大气压的比值；初始高度h₀；地面风速。

可选地，所述最大影响距离分别是在轻ERPGs的浓度值ERPG-3稳态下、中ERPGs的浓度值ERPG-2稳态下和重ERPGs的浓度值ERPG-1稳态下计算得到的。

可选地，所述本地坐标系是以下风向为x轴的正方向，垂直方向为y轴方向，泄漏源点的坐标设为(0，0)。

可选地，所述步骤S30中生成所述左边界包括：

S31，生成上风向左边界；

S32，生成下风向左边界以及末端左边界。

可选地，所述步骤S31的生成上风向左边界，具体包括：

步骤S311，从坐标(0，b₀)开始向上风向方向按照π/25的弧度间隔使用弧度推进法逆推计算5个点，每个所述点的坐标值为：

x＝pre_x-pre_y*sin(θ)； (8)

其中，所述点的x坐标值x表示上风向距离的相反数，y坐标值为b，表示重气云团半宽；g表示重力加速度，P₀表示重气密度，P_a表示空气密度；θ为弧度间隔值，pre_x表示前一点的x坐标值，pre_y表示前一点的y坐标值；

步骤S312，判断重气云团半宽b的大小是否大于5：如果大于5，则返回所述步骤S311；如果小于5，则向所述上风向方向按照π/10的弧度间隔使用半圆法逆推计算5个点，并按照所述公式(8)计算下风向距离x；所述公式(5)重气云团半宽b。

可选地，所述步骤S32的生成下风向左边界和末端左边界包括：

步骤S321，设定初始的步长step；

步骤S322，从所述本地坐标系的坐标(0，b₀)开始，沿着x轴的正方向，按照所述步长step递增取x值；并根据所述公式(5)计算重气云团半宽b；

步骤S323，判断x是否达到所述最大影响距离：如果达到，则将(x，b)作为所述下风向左边界的边界点，并跳转至所述步骤S325；如果没有达到，则跳转至步骤S324；

步骤S324，判断x是否达到风速*泄漏后的时间：如果达到，则跳转至步骤S328；如果未达到，则调整所述步长step，并跳转至所述步骤S322继续按照调整过的所述步长step进行递增取点；其中，调整的所述步长是依据下风向的边界距离而设定的；

步骤S325，使用所述弧度推进法从所述下风向左边界的边界点求取等值点坐标的x坐标值x，根据所述等值点的x坐标值x通过所述公式(5)计算重气云团半宽b，得到所述等值点坐标(x，b)；

步骤S326，判断所述等值点坐标的x坐标值x是否大于风速*泄漏后的时间：如果大于，则跳转至步骤S328；如果小于，则跳转至步骤S327；

步骤S327，判断所述重气云团半宽b是否收敛：如果收敛，则跳转至步骤S328；如果不收敛，则重新跳转回所述步骤S325；

步骤S328，采用所述半圆法求取末端等值线点坐标；

其中，所述末端等值线点为所述末端左边界的边界点。

可选地，所述步骤S40的校验和纠正所述左边界包括两部分：1)对所述ERPGs浓度为ERPG-3和ERPG-2时的两个等值线点的坐标集合进行校验和纠正；2)对所述ERPGs浓度为ERPG-2和ERPG-1的两个等值线点的坐标集合进行校验和纠正。

可选地，所述步骤S60的转换坐标系具体包括：

步骤S61，旋转所述左边界和所述右边界的边界点的本地坐标；

步骤S62，根据所述泄漏源的经纬度，计算出所述泄漏源的UTM坐标；

步骤S63，根据所述泄漏源的UTM坐标，将所述步骤S61旋转过的所述左边界和所述右边界的边界点的坐标转换为对应的UTM坐标；

步骤S64，将所述左边界和所述右边界的边界点的UTM坐标转换成经纬度坐标。

如上所述，本发明的一种重气持续泄漏扩散的动态模拟方法，具有以下有益效果：

1)优化了危害区域等浓度边界等值线，使得模拟所述的态势结果更接近实际的气体扩散过程；

2)优化了重气浓度的横向分布，使得重气浓度在横风向的分布上也有了一定的梯度变化；

3)本发明不仅仅可以模拟出稳定状态下的重气持续泄漏扩散态势，还可以动态的模拟重气持续泄漏扩散随时间变化的过程。

附图说明

图1显示为本发明的实施例公开的一种重气持续泄漏扩散的动态模拟方法的流程示意图。

图2显示为本发明的实施例公开的一种重气持续泄漏扩散的动态模拟方法中生成上风向左边界的流程示意图。

图3显示为本发明的实施例公开的一种重气持续泄漏扩散的动态模拟方法中生成下风向左边界和末端边界的流程示意图。

图4显示为本发明的实施例公开的一种重气持续泄漏扩散的动态模拟方法中对边界点进行坐标转换的流程示意图。

图5显示为在实施例2的已知参数下，使用平板模型模拟的ERPGs的等值线GIS效果图。

图6至图9显示为本发明实施例2公开的一种重气持续泄漏扩散的动态模拟方法模拟出的不同泄漏时间长度下的ERPGs的等值线GIS效果图。

元件标号说明

S10～S70 步骤

S31～S32 步骤

S311～S312 步骤

S321～S328 步骤

S61～S64 步骤

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

请参阅图1至图9，需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，遂图式中仅显示与本发明中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的型态、数量及比例可为一种随意的改变，且其组件布局型态也可能更为复杂。

平板模型是一种描述重气持续泄漏扩散的模型，具备较好的准确性和适度的计算复杂度，广泛用于突发事故的重气持续泄漏扩散的模拟；然而由于平板模型自身的一些数学上的局限性，采用一般的方法直接对模型进行求解模拟，效果比较差，对突发事故的指挥调度、救援等工作的指导意义不大。因此，本发明公开了一种基于平板模型的重气持续泄漏扩散的动态模拟方法。

实施例1

本实施例的一种重气持续泄漏扩散的动态模拟方法，具体如图1所示，包括：

步骤S10，输入已知参数，并根据已知参数初始化相关参数：

由于本实施例进行重气持续泄漏扩散的动态模拟的，所以其已知参数具体包括：10米高的风速u、风向(角度)、气压、气温、泄漏源经度long、泄露源纬度lat、泄漏源高度h、重气云团初始密度、重气云团初始半宽b₀、事发后的时间t、泄露物质分子量、泄漏物质的ERPGs浓度c和气云初始浓度c₀等等。

初始化相关参数具体包括：

a.将风向转化为弧度：弧度＝角度*π/180；

b.气压转化为与标准大气压的比值；

c.初始高度h₀：h₀＝b₀/2；

d.地面风速μ：其中，θ表示摩擦因子，在本实施例中，θ在0.2～0.4之间；

e.重气密度P₀：

步骤S20，计算ERPGs浓度稳态下的最大影响距离：

ERPGs是用于化工行业突发环境事件的评价标准，可以有效反映有毒化学品突发性泄漏对人体健康的影响程度。每种物质有分别代表危害程度轻、中、重的三种浓度值：ERPG-3、ERPG-2、ERPG-1。ERPGs浓度值是已知参数，带入平板模型：

V₀′＝2*b₀*h₀*u； (1)

h＝u_**u^-1B₀ ^-1/3*x^4/3； (3)

其中，V₀′表示初始气云体积通量；B₀表示浮力通量；g表示重力加速度，且g＝9.8m/s²；；P_a表示空气密度，且P_a＝1.2；u_*为摩擦风速，一般取10M高处的平均风速的6.67％；x为下风向距离；c₀表示气云初始浓度；c表示重气云浓度；b为重气云团半宽。

合并式(1)、(2)、(3)、(4)和(5)，可得出：

H*x³+x²＝temp； (6)

其中，H和temp是中间变量， 求解公式(5)，可以计算出最大影响距离S：

在本实施例中，根据ERPG-3、ERPG-2、ERPG-1三个轻、中、重三个浓度c，通过公式(7)计算得出三个浓度稳态下的最大影响距离S。

步骤S30，建立本地坐标系：以下风向为x轴的正方向，垂直风向为y轴方向，泄漏源点的坐标设为(0，0)；并在本地坐标系下生成左边界，具体包括：

步骤S31，生成上风向左边界；

步骤S32，生成下风向左边界和末端左边界。

其中，步骤S31中，生成上风向左边界具体如图2所示，包括：

步骤S311，使用弧度推进法从坐标(0，b₀)开始向上风向方向按照的π/25弧度间隔逆推计算5个点，即按照如下公式计算下风向距离x：

x＝pre_x-pre_y*sin(θ)； (8)

其中，θ为弧度间隔值，在本步骤中，由于使用的是弧度推进法，所以θ＝π/25；pre_x表示前一点的x坐标值；pre_y表示前一点的y坐标值；然后使用公式(5)计算重气云团半宽b，重气云团半宽b即纵坐标y；

步骤S312，判断重气云团半宽b的大小是否大于5：如果大于5，则返回步骤S311；如果小于5，则向上风向按照π/10弧度的间隔逆推计算5个点(半圆法)，即设置弧度间隔值θ＝π/10，按照式(8)计算下风向距离x，并使用公式(5)计算对应的重气云团半宽b；步骤S311和步骤S312中的所有点即为上风向左边界的所有边界点。

其中，步骤S32中生成下风向左边界和末端边界，具体如图3所示，包括：

步骤S321，设定初始的步长step为6M；

步骤S322，从坐标(0，b₀)开始，沿着x轴的正方向，即下风向方向，按照步长step递增取x值，即x＝x+step；带入公式(5)计算重气云团半宽b，从而得到一个等浓度值点；

步骤S323，判断此时的下风向距离是否达到步骤S20中的最大影响距离S：如果未达到，则跳转至步骤S324；如果达到，则将此时的点作为下风向左边界的边界点，并跳转至步骤S325；

步骤S324，判断此时的下风向距离是否达到u*t：如果达到，跳转至步骤S328；如果未达到，调整步长step，并跳转至步骤S322，继续按照调整过的步长step递增取点；其中，调整的步长step是依据下风向边界距离而设定的：如果下风向距离等于下风向边界距离，即120M、600M、1600M或4000M，对应的调整的步长step为：24M、50M、120M和300M；如果下风向距离不等于下风向边界距离120M、600M、1600M或4000M，则不对步长step进行调整；

步骤S325，使用弧度推进法从下风向左边界点开始，按照地推弧度θ＝π/25，求取等值点坐标的x坐标值x，根据等值点的x坐标值x通过公式(5)计算重气云团半宽b，得到等值点坐标(x，b)；

步骤S326，对于求取的所有等值点坐标，判断其下风向距离x是否大于u*t；如果大于，则跳转至步骤S328；如果小于，则跳转至步骤S327；

步骤S327，根据所有等值点的下风向距离x和公式(5)计算重气云团半宽b，并判断重气云团半宽距离b是否收敛：如果收敛，则跳转至步骤S328；如果不收敛，则重新跳转回步骤S325。

步骤S328，采用步骤S33中提到的半圆法求取末端等值线点坐标，并且，末端等值线点即为末端左边界的边界点。

步骤S40，为了防止发生边界的交叉和越界，本实施例还增加了对步骤S30生成的上风向左边界、下风向左边界和末端左边界的校验和纠正，包括：输入两个等值线点坐标的集合a和b，其中a是浓度为A的等值线点的坐标集合，b是浓度为B的等值线点的坐标集合，且A＜B。关于集合a和集合b的校验和纠正包括：提取集合a中的最大x坐标值max(x_a)，集合b中最大x坐标值max(x_b)，如果max(x_a)＜max(x_b)，采用补点的方法在集合a中增加新点，新点的x坐标的取值为最大坐标增加一定的步长，使得max(x_a)＞max(x_b)；如果max(x_a)＞max(x_b)，遍历集合a和集合b中的所有点坐标，判断集合a中是否存在任意点M，点M的x坐标值比集合b中的某点N的x坐标值大，但y坐标值却更小的点：如果存在，则调整M点的y坐标值为N点的y坐标值；如果不存在，则结束。

集合a包括{(x_a1,y_a1),(x_a2,y_a2),(...),(x_ai,y_ai),(...),(x_aI,y_aI)}，集合b包括{(x_b1,y_b1),(x_b2,y_b2),(...),(x_bj,y_bj),(...),(x_bJ,y_bJ)}。“遍历集合a和集合b中的所有点坐标，判断集合a中是否存在任意点M，点M的x坐标值比集合b中的某点N的x坐标值大，但y坐标值却更小的点：如果存在，则调整M点的y坐标值为N点的y坐标值；如果不存在，则结束”是按照如下进行的：

步骤一，初始化参数i＝1,j＝1；

步骤二，判断集合a中的i点和集合b中的j点是否满足条件x_ai＞x_bj，且y_ai＜y_bj：如果满足，则b集合中的j点的y坐标值替换a集合中i点的y坐标值，即y_ai＝y_bj，跳转至步骤三；如果不满足，则直接跳转至步骤三；

步骤三，判断i是否等于I：如果i＝I，则跳转至步骤四；如果i≠I，则重新跳转至步骤二；

步骤四，判断j是否等于J：如果j＝J，则跳转至步骤五；如果j≠J，则重新跳转至步骤二；

步骤五，结束。

在本实施例中，需要进行两次边界校验和纠正：1)对ERPGs浓度为ERPG-3和ERPG-2时的两个等值线点的坐标集合进行校验和纠正；2)对ERPGs浓度为ERPG-2和ERPG-1的两个等值线点的坐标集合进行校验和纠正。

步骤S50，生成右边界：由于平板模型是对称的，所以相较于左边界的所有边界点：右边界的边界点的x坐标的值与左边界的x坐标的值相同；右边界的边界点的y坐标的值与左边界的z坐标的值相反。

步骤S60，转换坐标系，具体如图4所示：

步骤S61，将上述步骤计算得到的边界点的坐标进行旋转：x轴正方向由风向转向正东方向：

X＝x*cos(n)+y*sin(n),Y＝y*cos(n)-x*sin(n)； (9)

步骤S62，根据由泄漏源经度long和泄漏源纬度lat组成的泄漏源经纬度坐标O(long，lat)，按照如下公式计算出泄漏源经纬度坐标对应的UTM坐标O(x，y，xzone，yzone)：

1)计算xzone：如果long＜0，那么如果long≥0，那么

2)计算yzone：查找lat值在数组【-90,-84,-72,-64,-56,-48,-40,-32,-24,-16,-8,0,8,16,24,32,40,48,56,64,72,84】中的位置，与数组【'A','C','D','E','F','G','H','J','K','L','M','N','P','Q','R','S','T','U','V','W','X','Z'】相对应位置的值即为yzone的值；

3)计算x和y：

n＝(a-b)/(a+b)；

rh₀＝a(1-e²)/(1-e²sin²(lat))^3/2；

n_u＝a/(1-e²sin²(lat))^1/2；

long₀＝6*xzone-183；

p＝(long-long₀)；

M＝a[(1-e²/4-3e⁴/64-5e⁶/256)lat

-(3e²/8+3e⁴/32+45e⁶/1024)sin(2lat)

+(15e⁴/256+45e⁶/1024)sin(4lat)

-(35e⁶/3072)sin(6lat))]；

y＝K₁+K₂p²+K₃p⁴；

K₁＝M*K₀；

K₂＝K₀n_usin(lat)cos(lat)/2＝K₀n_usin(2lat)/4；

K₃＝[K₀n_usin(lat)cos³(lat)/24][5-tan²(lat)+

9*elsq*cos²(lat)+4*elsq²*cos⁴(lat)]；

x＝K₄+K₅p³；

K₄＝K₀n_ucos(lat)；

K₅＝(K₀n_ucos³(lat)/6)[1-tan²(lat)+

elsq*cos²(lat)]。

其中，a、b、e、K₀和elsq均为常数，且a＝6378137；b＝6356752.314；e＝0.081819191；K₀＝0.9996；elsq＝0.006739497；

步骤S63，结合泄漏源点在本地坐标系中的坐标(0，b₀)，可分别求出经过步骤S61旋转后的边界点的坐标对应的UTM坐标：

源点的UTM坐标为(x0，y0，xzone，yzone)，那么Z(x，y)的UTM坐标就为：(x0+x，y0+y，xzone，yzone)；

判断是否发生越区，即x₀+x＜0或x₀+x＞1000000：如果发现越区，则xzone±1；x₀+x±500000；

步骤S64，将边界点的UTM坐标转换成经纬度坐标：

根据等值线点UTM坐标的yzone值来判定等值线点位于哪个半球：yzone在纬度区A、C、D、E、F、G、H、J、K、L和M位于南半球，而其余区位于北半球。

按照如下公式计算经纬度坐标：

arc＝y/k₀；

m_u＝arc/[a*(1-e²/4-3*e⁴/64-5*e⁶/256)]；

e_i＝[1-(1-e²)^1/2]/[1+(1-e²)^1/2]；

c_a＝(3*e_i/2-27*e_i ³/32)；

c_b＝(21*e_i ²/16-55*e_i ⁴/32)；

c_c＝(151*e_i ³/96)；

c_d＝(1097*e_i ⁴/512)；

phi1＝m_u+c_a*sin(2m_u)+c_b*sin(4*m_u)+

c_c*sin(6*m_u)+c_d*sin(8*m_u)；

Q₀＝elsq*cos²(phi1)；

t₀＝tan²(phi1)；

r₀＝a(1-e²)/(1-e²sin²(phi1))^3/2；

n₀＝a/(1-e²sin²(phi1))^1/2；

dd₀＝x/(n₀*k₀)；

long₀＝6*xzone-183；

fact₁＝n₀*tan(phi1)/r₀；

fact₂＝(dd₀)²/2；

fact₃＝(5+3*t₀+10*Q₀-4*Q₀ ²-9*elsq)*(dd₀)⁴/24；

fact₄＝(61+90t₀+298*Q₀+45*(t₀)²-3*(Q₀)²-252*elsq)*(dd₀)⁶/720；

lat＝[phi1-fact₁(fact₂+fact₃+fact₄)]；(如果是南半球，则lat＝-lat)；

loft₁＝dd₀；

loft₂＝(1+2*t₀+Q₀)*(dd₀)³/6；

loft₃＝(5-2*Q₀+28*t₀-3*(Q₀)²+8*elsq+

24*(t₀)²)*(dd₀)⁵/120；

long＝long₀-(loft₁-loft₂+loft₃)*cos(phi1)。

步骤S70，输出所有的等值线点的经纬度坐标，并可根据该经纬度坐标绘制GIS(Geographic Information System或Geo－Information system，地理信息系统)图层，显示于GIS平台上。

实施例2

为了更清楚地说明本发明，本实施例将初始参数设置为：

在初始参数为上表所述的情况下，直接应用平板模型求解得到的ERPGs的等值线GIS效果图如图5所示。使用本实施例的方法，求解的ERPGs的等值线效果如下：

1)泄漏源泄漏5分钟后，ERPGs的等值线GIS效果图如图6所示；

2)泄漏源泄漏10分钟后，ERPGs的等值线GIS效果图如图7所示；

3)泄漏源泄漏20分钟后，ERPGs的等值线GIS效果图如图8所示；

4)泄漏源泄漏30分钟后，ERPGs的等值线GIS效果图如图9所示。

从图5～图9不难看出，使用本实施例的方法，可得到更好的动态模拟结果：优化了危害区域边界等值线，使得模拟所得的态势结果更接近实际的重气扩散过程；优化了重气浓度在横风向上的分布，使得其在横风向上有了一定的梯度变化。

此外，为了使本实施例所公开的方法能够更方便的被第三方使用，已将该方法基于java语言实现，并封装成了jar包：sari.sel.model.jar，其具体实施步骤如下：

1，导入jar包；将sari.sel.model.jar导入到目标java工程目录下；

2，在需要引用该方法的类文件的头部，加入import sari.sel.model.PingBanModel，import sari.sel.model.CoordinatePoints.GPS_Point；

3，模型类的实例化；PingBanModel mod＝new PingBanModel()；

4，通过set方法，设置实例的具体参数(包括：风速、风向、气温、气压、初始浓度、重气云团初始半宽、泄漏源经纬度、泄漏源高度、泄漏物质分子量、泄漏物质对应的ERPGs浓度值，仿真时间)；设置方法如下，例如：mod.setWind_Speed(5),设置风速大小为5m/s；

5，调用mod.Start_Simulate()方法，开始仿真计算；

6，分别调用mod.getGPS_Coordinates_Low()，mod.getGPS_Coordinates_Low()，mod.getGPS_Coordinates_Low()三个方法，可分别得到该泄露物质ERPGs的三个浓度，分别对应的等浓度线上的点经纬度坐标集合，返回结果的数据类型为ArrayList<GPS_Point>，其中GPS_Point为自定义类，包括lon，lat两个属性。

综上所述，本发明的一种重气持续泄漏扩散的动态模拟方法，具有以下有益效果：优化了末端等值线，使得模拟所得的态势结果更接近实际的气体扩散过程；优化了重气浓度的横向分布，使得重气在横向上也有了一定的梯度变化；本发明不仅仅可以模拟出稳定状态下的重气持续泄漏扩散态势，还可以动态的模拟重气持续泄漏扩散随时间变化的过程。所以，本发明有效克服了现有技术中的种种缺点而具高度产业利用价值。

上述实施例仅例示性说明本发明的原理及其功效，而非用于限制本发明。任何熟悉此技术的人士皆可在不违背本发明的精神及范畴下，对上述实施例进行修饰或改变。因此，举凡所属技术领域中具有通常知识者在未脱离本发明所揭示的精神与技术思想下所完成的一切等效修饰或改变，仍应由本发明的权利要求所涵盖。

Claims (7)

1.一种重气持续泄漏扩散的动态模拟方法，所述方法是基于平板模型的，其特征在于，所述动态模拟方法包括：

步骤S10，输入已知参数，并根据所述已知参数初始化相关参数；

步骤S20，计算ERPGs浓度稳态下的最大影响距离；

步骤S30，在本地坐标系下生成左边界：其中，所述本地坐标系是以下风向为x轴的正方向，垂直方向为y轴方向，泄漏源点的坐标设为(0，0)；包括：S31，生成上风向左边界；S32，生成下风向左边界以及末端左边界；所述步骤S31具体包括：步骤S311，从坐标(0，b₀)开始向上风向方向按照π/25的弧度间隔使用弧度推进法逆推计算5个点，每个所述点的坐标值为：

$b=b_0*{(1+1.5*{\left(\frac{g*h_0*(P_0-P_a)}{P_a}\right)}^{\frac{1}{2}}*\frac{x}{u*b_0})}^{\frac{2}{3}};---\left(5\right)$

x＝pre_x-pre_y*sin(θ)； (8)

其中，所述点的x坐标值x表示上风向距离的相反数，y坐标值为b，表示重气云团半宽；g表示重力加速度，P_a表示空气密度；θ为弧度间隔值，pre_x表示前一点的x坐标值，pre_y表示前一点的y坐标值；步骤S312，判断重气云团半宽b的大小是否大于5：如果大于5，则返回所述步骤S311；如果小于5，则向所述上风向方向按照π/10的弧度间隔使用半圆法逆推计算5个点，并按照所述公式(8)计算下风向距离x；所述公式(5)计算重气云团半宽b；

步骤S40，校验纠正所述左边界；

步骤S50，生成右边界：右边界的边界点的x坐标的值与所述左边界的边界点的x坐标的值相同，y坐标的值相反；

步骤S60，转换坐标系：将所述左边界和所述右边界的边界点的坐标转换为经纬度坐标；

步骤S70，输出所述步骤S60转换得到的经纬度坐标，并据此绘制GIS效果图。

2.根据权利要求1所述的重气持续泄漏扩散的动态模拟方法，其特征在于，所述已知参数包括风速u、风向、气压、气温、泄漏源的经纬度、泄漏源的高度、重气云团初始密度、重气云团初始半宽b₀、泄漏后的时间、重气的分子量和ERPGs浓度。

3.根据权利要求2所述的重气持续泄漏扩散的动态模拟方法，其特征在于，所述相关参数包括：用弧度表示的风向；气压与标准大气压的比值；初始高度h₀；地面风速和重气密度P₀。

4.根据权利要求3所述的重气持续泄漏扩散的动态模拟方法，其特征在于，所述最大影响距离分别是在轻ERPGs的浓度值ERPG-3稳态下、中ERPGs的浓度值ERPG-2稳态下和重ERPGs的浓度值ERPG-1稳态下计算得到的。

5.根据权利要求1所述的重气持续泄漏扩散的动态模拟方法，其特征在于，所述步骤S32的生成下风向左边界和末端左边界包括：

步骤S321，设定初始的步长step；

步骤S322，从所述本地坐标系的坐标(0，b₀)开始，沿着x轴的正方向，按照所述步长step递增取x值；并根据所述公式(5)计算重气云团半宽b；

步骤S323，判断x是否达到所述最大影响距离：如果达到，则将(x，b)作为所述下风向左边界的边界点，并跳转至步骤S325；如果没有达到，则跳转至步骤S324；

步骤S324，判断x是否达到风速*泄漏后的时间：如果达到，则跳转至步骤S328；如果未达到，则调整所述步长step，并跳转至所述步骤S322继续按照调整过的所述步长step进行递增取点；其中，调整的所述步长是依据下风向的边界距离而设定的；

步骤S325，使用所述弧度推进法从所述下风向左边界的边界点求取等值点坐标的x坐标值x，根据所述等值点的x坐标值x通过所述公式(5)计算重气云团半宽b，得到所述等值点坐标(x，b)；

步骤S326，判断所述等值点坐标的x坐标值x是否大于风速*泄漏后的时间：如果大于，则跳转至步骤S328；如果小于，则跳转至步骤S327；

步骤S327，判断所述重气云团半宽b是否收敛：如果收敛，则跳转至步骤S328；如果不收敛，则重新跳转回所述步骤S325；

步骤S328，采用所述半圆法求取末端等值线点坐标；

其中，所述末端等值线点为所述末端左边界的边界点。

6.根据权利要求4所述的重气持续泄漏扩散的动态模拟方法，其特征在于，所述步骤S40的校验和纠正所述左边界包括两部分：1)对所述ERPGs浓度为ERPG-3和ERPG-2时的两个等值线点的坐标集合进行校验和纠正；2)对所述ERPGs浓度为ERPG-2和ERPG-1的两个等值线点的坐标集合进行校验和纠正。

7.根据权利要求6所述的重气持续泄漏扩散的动态模拟方法，其特征在于：所述步骤S60的转换坐标系具体包括：

步骤S61，旋转所述左边界和所述右边界的边界点的本地坐标；

步骤S62，根据所述泄漏源的经纬度，计算出所述泄漏源的UTM坐标；

步骤S63，根据所述泄漏源的UTM坐标，将所述步骤S61旋转过的所述左边界和所述右边界的边界点的坐标转换为对应的UTM坐标；

步骤S64，将所述左边界和所述右边界的边界点的UTM坐标转换成经纬度坐标。