CN104657573A - 用于三维空间的泄漏气体扩散预测方法 - Google Patents

Abstract

本发明提供了一种用于三维空间的泄漏气体扩散预测方法。该方法通过构建空间三维坐标系，将泄漏区域及周边区域表示在同一坐标系中，并把其地形及建筑物分布数据加载在该坐标系中的固定区计算域内进行泄漏气体扩散模拟，在此基础上进行固定区计算域外的毒气扩散预测，从而实现泄漏源周边地区的泄漏气体扩散浓度分布预测。本发明在对泄漏气体进行扩散预测时充分考虑当时当地的地理信息，并能够快速计算出一个或多个泄漏源造成的污染区域的范围和被污染区域内的毒气浓度。

Description

用于三维空间的泄漏气体扩散预测方法

技术领域

本发明涉及气体扩散预测领域。

背景技术

2006年2月，国务院发布《国家中长期科学和技术发展规划纲要（2006～2020年）》，将“公共安全”列为重点领域，明确提出发展危险化学品泄漏应急救援技术、核生化恐怖源的远程探测技术等。2011年3月，时任国务院总理温家宝签署国务院令公布修订后的《危险化学品安全管理条例》。条例中所涉及的危险化学品，特别是有毒气体的泄漏，会在风力的作用下可快速传播至几公里之外，导致大面积中毒，对人体、设施、环境等都会造成危害，甚至是永久性伤害。因此，在控制泄漏源的同时，需要迅速疏散泄漏源周围的居民，并进行中毒人员救援，而这些都需要首先确定泄漏气体的扩散范围、趋势，从而确定警戒与疏散范围、撤离路线以及救援资源分配等。上述这些工作均需要建立在对泄漏气体的实时扩散进行模拟的基础上，如何合理地、有效地通过数值模拟预测泄漏物质的迁移演化规律有利于社会安全，也是提升应急救援效率的重要前提。

目前，常用的扩散模型：如高斯模型（Gaussian）、高斯烟羽模型（Gaussianplume）、高斯烟团模型（Gaussian puff）、BM模型、Sutton模型、三位有限元计算模型（FEM3）等，这些模型一定程度为泄漏气体扩散预测奠定了基础，国内学者丁信伟等（丁信伟，王淑兰.可燃及毒性气体泄漏扩散研究综述.化学工业与工程，1999，16(2):118-122）对以上几类模型进行分析并指出其各自的适用范围和局限，在此不再叙述。在上述所有模型中，高斯模型提出时间较早、实验数据多并且计算方便。目前，很多应用都采用高斯模型或其改进模型，但由于实际环境较为复杂，特别是地质环境（如实际扩散区域大多有起伏地形、各类建筑等多种地形因素）以及泄漏源数量等严重影响了实际预测效果。

2009年，Fan等（Fan L,Yi S,Li Z,et al.Research on toxic Gas DiffusionSimulation in Urban areas based on GIS.New Trends in Information and ServiceScience,NISS′09.International Conference on.IEEE,2009:166-169）利用改进的高斯模型模拟城区内的气体扩散，并达到一定精度，但其对于扩散预测效率和输入条件方面的要求较为严格，实际应用中会受到一定限制；同时，也未考虑到同时有2个或2个以上泄漏源时所带来的影响。2011年，Stockie（Stockie J M.The mathematics of atmospheric dispersion modeling.SiamReview,2011,53(2):349-372）重现了高斯模型的推导过程，根据物理假设的不同，得出不同的结果，但仍旧未考虑空间三维地形以及泄漏源数量等的影响。总的来说，对具有复杂几何学特征的高密度建筑物区域以及高起伏地形条件下的多泄漏源扩散浓度值的预测，高斯模型的性能并不理想。

2011年，黄弘等提出了一种基于Lagrangian模型与Eulerian模型耦合的建筑物周边气体扩散模拟方法（黄弘，胡啸峰等，基于Lagrangian模型与Eulerian模型耦合的建筑物周边气体扩散模拟,清华大学学报（自然科学版），2011，51(12):1870-1876），但该方法未涉及高起伏地形条件以及多泄漏源扩散问题，该方法主要适用于城市建筑物周边大气扩散预测。

随着计算机的计算性能提升，基于计算流体力学模型的方法得到了很大发展（Zajaczkowski F J,Haupt S E,Schmehl K J.A preliminary study ofassimilating numerical weather prediction data into computational fluid dynamicsmodels for wind prediction.Journal of Wind Engineering and IndustrialAerodynamics,2011,99(4):320-329），并成为近年来人们研究泄漏气体扩散的主要趋势，但由于计算流体力学模型需要对平面和空间进行离散化，尽管计算性能提升，但要实现大范围，特别是高起伏地形有毒气体扩散预测（相对于平面和平坦三维地形，高起伏地形条件会引起计算域空间突增），难度十分大，其主要原因在于对整个关注区域都建模计算，这无疑增大了求解的难度，降低了扩散求解的效率，不利于实际应用。

发明内容

本发明提供了一种用于三维空间的泄漏气体扩散预测方法，该方法能对具有复杂几何学特征的高密度建筑物区域和/或高起伏地形下的一个或多个泄漏源形成的扩散浓度场进行预测，并兼顾预测效率。

根据本发明的一个方面，提供了一种用于三维空间的泄漏气体扩散预测方法，该方法包括：步骤S1，采集参数，所述参数包括地理环境参数和泄漏源自身的参数；步骤S2，建立坐标系，并参照该坐标系建立固定区计算域，将所述固定区计算域内的气体离散化为由多个离散点界定的流体微团；步骤S3，基于流体力学领域的纳维-斯托克斯方程建立在所述固定区计算域内的泄漏气体扩散方程，所述泄漏气体扩散方程中包括质量守恒方程、动量方程和能量守恒方程；步骤S4，将所述步骤S3中得到的所述质量守恒方程、所述动量方程和所述能量守恒方程转换为关于空间离散坐标和时间坐标的差分方程，并根据定解条件和采集到的所述参数，以迭代的方式求解所述固定区计算域内的泄漏气体扩散浓度分布场的数值解；步骤S5，建立所述泄漏气体扩散浓度分布场与采集到的所述参数中的部分或全部间的函数关系；步骤S6，将所述步骤S5中得到的所述函数关系表示为包含待定系数的无量纲化函数关系，并且将从所述步骤S4中得到的所述数值解中的部分或全部代入所述无量纲化函数关系，以确定所述待定系数的值；以及步骤S7，将感兴趣的点的空间坐标和时间坐标代入步骤S6中得到的所述系数已被确定的所述无量纲化函数关系，计算得到所述感兴趣的点的泄漏气体扩散浓度，所述感兴趣的点在所述固定区计算域内或者在所述固定区计算域外。

当有多个泄漏源时，针对所述多个泄漏源中的每一个，按照所述步骤1～步骤7计算该泄漏源在所述感兴趣的点的所述泄漏气体扩散浓度，将计算得到的所述多个泄漏源的所述泄漏气体扩散浓度相加，以得到所述感兴趣的点的总的泄漏气体扩散浓度；或者针对所述多个泄漏源中的每一个，按照所述步骤1～步骤6得到对应于该泄漏源的所述泄漏气体扩散浓度分布场的所述无量纲化函数关系，总的泄漏气体扩散浓度分布场是所述多个泄漏源的所述泄漏气体扩散浓度分布场的和，再将所述感兴趣的点的空间坐标和时间坐标代入所述总的泄漏气体扩散浓度分布场的所述无量纲化函数关系，以得到所述感兴趣的点的总的泄漏气体扩散浓度。

本发明通过构建空间三维坐标系，将毒气泄漏区域及周边区域表示在同一坐标系中，并把其地形及建筑物分布数据加载在坐标系中的一定区域内进行毒气扩散模拟，在此基础上进行固定区域外的毒气扩散预测，从而实现泄漏源周边地区的毒气浓度分布预测。本发明所提供的方法能在对泄漏气体进行扩散预测时充分考虑当时当地的地理信息，并能够快速计算出一个或多个泄漏源造成的污染区域的范围和被污染区域内的毒气浓度。

附图说明

图1示出了示出了根据本发明的优选实施方式的一种用于三维空间的泄漏气体扩散预测方法；以及

图2示出了图1所示的方法中对固定区计算域进行离散化的示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的附图标记。附图中未绘示或描述的实现方法，为所属技术领域中普通技术人员所知的方法。另外，虽然本文可提供包含特定值的参数的示范，但应该了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。此外，以下实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向。因此，使用的方向用语是用来说明并非用来限制本发明。

图1示出了示出了根据本发明的优选实施方式的一种用于三维空间的泄漏气体扩散预测方法。

步骤S1：采集所需的参数，采集的参数中可包括地理环境参数和泄漏源自身的参数，地理环境参数可包括泄漏源所在地的经度、纬度和海拔(X₀,Y₀,H₀)（经度和纬度的单位可为°，海拔的单位可为m）、云层密度ρ₁（Kg/m³）、风速V（m/s）、风向α（为风向与正东向之间的夹角，单位为弧度）和温度T（℃）等，泄漏源自身的参数包括泄漏气体密度ρ₀（ug/m³）、泄漏类型参数（若为瞬时泄漏为泄漏量Q（Kg/s）、若为连续泄漏则为泄漏时间T_t0（s）和泄漏质量流率q（Kg/s））和泄漏源特征尺寸L（m）（泄漏源特征尺寸为泄漏源的三维结构沿某一方向上的最大尺寸值）。

泄漏源经度、纬度和海拔(X₀,Y₀,H₀)可以通过GPS（全球定位系统）接收机进行测量，也可以通过手机定位和地图进行数值提取。

风速V，温度T可以通过气象数据提取，也可以通过现场测量获得。

还可根据需要采集泄漏气体扩散控制参数K(X₀,Y₀,H₀)等，泄漏气体扩散控制参数K(X₀,Y₀,H₀)可根据经验进行选择，通常选择范围为0至1间的某个数。

泄漏源特征尺寸L可通过激光测距仪获取，也可以通过目测进行估计。

可使用风向标对风向进行精确测量，或从天气预报参数中提取风向。可设α∈[0,2π)，用来表示风向在0至360度范围内取值。可对风向α进行如下划分：

步骤S2：建立坐标系，并参照该坐标系建立固定区计算域，将所述固定区计算域内的气体离散化为多个离散点界定的流体微团。

步骤S21：如图2所示，为使方程尽量简洁，优选地，可以以泄漏源中心为坐标原点O，以泄漏源大气风向作为直角坐标系的X正方向，建立三维右手坐标系OXYZ。

步骤S22：如图2所示，以泄漏源特征尺寸L为基准，建立以坐标原点O为中心的正方体计算域，坐标原点O到正方体各平面之间的垂直距离为βL，其中，β应为远大于1，且小于等于1000的一个数值，相应的计算域为(2βL)×(2βL)×(2βL)的正方体。

步骤S23：如图2所示，对步骤S22中得到的正方体计算域沿X轴、Y轴和Z轴进行离散化，离散化后的间隔分别为Δx，Δy和Δz，相应的X轴、Y轴和Z轴分别对应的点数为 和其中，表示取小于等于数值^*的最大整数，相应的计算域空间坐标为(x,y,z)，其中 $x=\mathrm{\Δx}\×(-\frac{N_x}{2}:1:(\frac{N_x}{2}-1)),$ $y=\mathrm{\Δy}\×(-\frac{N_y}{2}:1:(\frac{N_y}{2}-1)),$ $z=\mathrm{\Δz}\×(-\frac{N_z}{2}:1:(\frac{N_z}{2}-1)).$ 每个体积为Δx×Δy×Δz的空间内的气体可被视为一个流体微团。

实际区域中的泄漏气体通常存在于地表上方，而极少甚至不会扩散到地表下方。实际离散过程中，z的取值能够体现地形起伏或建筑物等因素。因此，本方法可以很好地应用于具有复杂几何学特征的高密度建筑物区域和高起伏地形条件下的扩散浓度值预测。

步骤S3：基于流体力学领域的纳维-斯托克斯（N-S）方程建立所述固定区计算域内的泄漏气体扩散方程。

由于泄漏气体在空气中的扩散是一种连续的流体介质运动过程，因此，其遵循流体力学的基本方程——N-S（纳维-斯托克斯）方程。本发明权衡计算精度和计算复杂度，对所研究的问题进行适当简化，来建立流体力学纳维-斯托克斯方程。

设固定区计算域内有N_Atm种组份（即组成成分，如空气和泄漏气体），其中将空气视为一种组分，而不考虑空气本身的构成。对于系统中的任一流体微团，其在运动过程中满足质量守恒，因此对于任意一种组份i_Atm，均满足式（1）：

$\frac{{\∂\mathrm{\ρM}}_{i_{\mathrm{Atm}}}}{\∂t}+\▿\·\left({\mathrm{\ρM}}_{i_{\mathrm{Atm}}}V\right)+\▿\·J_{i_{\mathrm{Atm}}}=R_{i_{\mathrm{Atm}}}+S_{i_{\mathrm{Atm}}}---\left(1\right)$

其中，▽为拉普拉斯算子，ρ（Kg/m³）为当地混合气体密度（可取固定区计算域内的混合气体密度的平均值或其中某一点的气体密度）；为组份i_Atm的浓度（即质量百分比，组份i_Atm的质量除以混合气体的质量的百分比）；V为所取流体微团的运动速度矢量（m/s）；与分别为化学反应导致的组份i_Atm的速率和源项导致的组份i_Atm的速率，单位为Kg/(s·m³)；为组份i_Atm的扩散通量，由浓度梯度产生，在此问题中可以将扩散视为稀释，则有J_i＝-ρD_i,m▽M_i，其中D_i,m为组份i_Atm的扩散系数，通常可从气体的物质属性表中查得。

式（1）实际上应为一组方程，当系统中有N_Atm种组分时，组中方程的个数应为N_Atm-1，显然，当只有一种泄漏气体时，只需要建一个形式如式（1）的方程，且为了保证计算精度应先求解（此处i_Atm=2用于表示空气），最后由得到泄漏气体的浓度。

对于大范围开放空间内的泄漏气体扩散问题，需要考虑气体的彻体力，即需要考虑空气中的浮力效应。对于所关注区域内的任一流体微团，其在运动过程中满足动量方程，则有：

$\frac{\∂\left(\mathrm{\ρu}\right)}{\∂t}+\▿\·\left(\mathrm{\ρuV}\right)=-\frac{\∂p}{\∂x}+\frac{{\∂\mathrm{\τ}}_{\mathrm{xx}}}{\∂x}+\frac{{\∂\mathrm{\τ}}_{\mathrm{yx}}}{\∂y}+\frac{{\∂\mathrm{\τ}}_{\mathrm{zx}}}{\∂z}+{\mathrm{\ρf}}_x---\left(2\right)$

$\frac{\∂\left(\mathrm{\ρv}\right)}{\∂t}+\▿\·\left(\mathrm{\ρvV}\right)=-\frac{\∂p}{\∂y}+\frac{{\∂\mathrm{\τ}}_{\mathrm{xy}}}{\∂x}+\frac{{\∂\mathrm{\τ}}_{\mathrm{yy}}}{\∂y}+\frac{{\∂\mathrm{\τ}}_{\mathrm{zy}}}{\∂z}+{\mathrm{\ρf}}_y---\left(3\right)$

$\frac{\∂\left(\mathrm{\ρw}\right)}{\∂t}+\▿\·\left(\mathrm{\ρwV}\right)=-\frac{\∂p}{\∂z}+\frac{{\∂\mathrm{\τ}}_{\mathrm{xz}}}{\∂x}+\frac{{\∂\mathrm{\τ}}_{\mathrm{yz}}}{\∂y}+\frac{{\∂\mathrm{\τ}}_{\mathrm{zz}}}{\∂z}+{\mathrm{\ρf}}_z---\left(4\right)$

由于扩散带来的掺混，实际中气体密度会随时间和空间位置发生变化的，所以为了求得空间中的气体密度，控制方程中需要包含能量方程。假设泄漏气体在扩散过程中不包含焓的输运，故能量方程中不需要考虑物质输送。流体微团在运动过程中满足能量守恒方程：

$\frac{\∂}{\∂t}\left[\mathrm{\ρ}(e+\frac{V^2}{2})\right]+\▿\·\left[\mathrm{\ρ}(e+\frac{V^2}{2})V\right]$

$=\mathrm{\ρq}+\frac{\∂}{\∂x}\left(k\frac{\∂T}{\∂x}\right)+\frac{\∂}{\∂y}\left(k\frac{\∂T}{\∂y}\right)+\frac{\∂}{\∂z}\left(k\frac{\∂T}{\∂z}\right)-\frac{\∂\left(\mathrm{up}\right)}{\∂x}-\frac{\∂\left(\mathrm{vp}\right)}{\∂y}-\frac{\∂\left(\mathrm{wp}\right)}{\∂z}+$

$\frac{\∂\left({\mathrm{u\τ}}_{\mathrm{xx}}\right)}{\∂x}+\frac{\∂\left({\mathrm{u\τ}}_{\mathrm{yx}}\right)}{\∂y}+\frac{\∂\left({\mathrm{u\τ}}_{\mathrm{zx}}\right)}{\∂z}+\frac{\∂\left({\mathrm{v\τ}}_{\mathrm{xy}}\right)}{\∂x}+\frac{\∂\left({\mathrm{v\τ}}_{\mathrm{yy}}\right)}{\∂y}+\frac{\∂\left({\mathrm{v\τ}}_{\mathrm{zy}}\right)}{\∂z}+$

$\frac{\∂\left({\mathrm{w\τ}}_{\mathrm{xz}}\right)}{\∂x}+\frac{\∂\left({\mathrm{w\τ}}_{\mathrm{yz}}\right)}{\∂y}+\frac{\∂\left({\mathrm{w\τ}}_{\mathrm{zz}}\right)}{\∂z}+\mathrm{\ρf}\·V---\left(5\right)$

其中，u、v、w分别为速度矢量V在坐标轴的X、Y和Z三个方向的分量，单位为m/s；p为当地压强，单位为Pa，通常可以从气象信息中获得；f_x、f_y和f_z分别表示为作用在流体微团上的彻体力，宏观上体现为浮力；τ_xx、τ_yy、τ_zz、τ_xy、τ_yx、τ_xz、τ_zx、τ_yz和τ_zy为作用在流体微团上的各个方向上的粘性力（分子粘性力）；e指气体内能，其定义可参见下式（7），单位为焦耳；k是绝热指数，此处可取1.4。彻体力和粘性力的值在迭代计算的过程中得到。需要说明的是，式（1）～式（5）为一般流体流动中所遵循的方程，在本发明的这一优选实施方式中，为了平衡计算量和计算结果的精度，在求解前对所研究问题进行了合理简化，将式（1）～式（5）简化为下面式（8）的形式。

本发明中，空间中各处均不存在气体源项，所以可以对上述质量守恒方程（1）进行简化。可对固定区计算域内的气体做层流假设，由此可以对式（2～5）进行简化，将其转化为欧拉方程的形式。还可以假设流体微团在运动过程中是绝能的，即不存在热交换。

另外气体在运动过程中具有如式（6）的本构方程：

p＝ρRT   （6）

式中，ρ（Kg/m³）为当地混合气体密度，p（Pa）为当地压强，R为气体常数，T为温度，当温度不变时，从式（6）中可以得到p和ρ的关系。

根据上述方式和所做的简化假设，得到方程组（7）和方程组（8）：

$J_{i_{\mathrm{Atm}}}=-D_{i_{\mathrm{Atm}},m}\▿M_{i_{\mathrm{Atm}}};$

p＝ρRT;   （7）

e＝c_vT;

$\left\{\begin{array}{c}\frac{\∂\left({\mathrm{\ρM}}_{i_{\mathrm{Atm}}}\right)}{\∂t}+\▿\·\left({\mathrm{\ρM}}_{i_{\mathrm{Atm}}}V\right)+\▿\·J_{i_{\mathrm{Atm}}}=0;\\ \frac{\∂\left(\mathrm{\ρu}\right)}{\∂t}+\▿\·\left(\mathrm{\ρuV}\right)=-\frac{\∂p}{\∂x}+{\mathrm{\ρf}}_x;\\ \frac{\∂\left(\mathrm{\ρv}\right)}{\∂t}+\▿\·\left(\mathrm{\ρvV}\right)=-\frac{\∂p}{\∂y}+{\mathrm{\ρf}}_y;\\ \frac{\∂\left(\mathrm{\ρw}\right)}{\∂t}+\▿\·\left(\mathrm{\ρwV}\right)=-\frac{\∂p}{\∂z}+{\mathrm{\ρf}}_z;\\ \frac{\∂}{\∂t}\left[\mathrm{\ρ}(e+\frac{V^2}{2})\right]+\▿\·\left[\mathrm{\ρ}(e+\frac{V^2}{2})V\right]=-\frac{\∂\left(\mathrm{up}\right)}{\∂x}-\frac{\∂\left(\mathrm{vp}\right)}{\∂y}-\frac{\∂\left(\mathrm{wp}\right)}{\∂z}+\mathrm{\ρf}\·V\end{array}\right.---\left(8\right)$

为了方便表述，可将简化后的扩散控制方程组（8）写成向量形式：

$\frac{\∂U}{\∂t}=J-\frac{\∂F}{\∂x}-\frac{\∂G}{\∂y}-\frac{\∂H}{\∂z}---\left(9\right)$

其中， $U=\left\{\begin{array}{c}{\mathrm{\ρM}}_{i_{\mathrm{Atm}}}\\ \mathrm{\ρu}\\ \mathrm{\ρv}\\ \mathrm{\ρw}\\ \mathrm{\ρ}(e+\frac{V^2}{2})\end{array}\right\},$ $F=\left\{\begin{array}{c}{\mathrm{\ρM}}_{i_{\mathrm{Atm}}}u+J_{i_{\mathrm{Am}}}^x\\ {\mathrm{\ρu}}^2+p\\ \mathrm{\ρuv}\\ \mathrm{\ρuw}\\ \mathrm{\ρu}(e+\frac{V^2}{2})+\mathrm{pu}\end{array}\right\},$ $F=\left\{\begin{array}{c}{\mathrm{\ρM}}_{i_{\mathrm{Atm}}}v+J_{i_{\mathrm{Am}}}^y\\ \mathrm{\ρvu}\\ {\mathrm{\ρv}}^2+p\\ \mathrm{\ρvw}\\ \mathrm{\ρv}(e+\frac{V^2}{2})+\mathrm{pv}\end{array}\right\},$ $H=\left\{\begin{array}{c}{\mathrm{\ρM}}_{i_{\mathrm{Atm}}}w+J_{i_{\mathrm{Am}}}^z\\ \mathrm{\ρwu}\\ \mathrm{\ρwv}\\ {\mathrm{\ρw}}^2+p\\ \mathrm{\ρw}(e+\frac{V^2}{2})+\mathrm{pw}\end{array}\right\},$ $J=\left\{\begin{array}{c}0\\ {\mathrm{\ρf}}_x\\ \mathrm{\ρ}{f}_y\\ {\mathrm{\ρf}}_z\\ \mathrm{\ρ}({\mathrm{uf}}_x+{\mathrm{vf}}_y+{\mathrm{wf}}_z)\end{array}\right\}.$

方程（7）和（8）或者方程（7）和（9）可被认为是泄漏气体扩散方程。

步骤S4：求解固定区计算域内的泄漏气体扩散浓度分布场。上述方程组为非线性偏微分方程组，目前尚不存在解析形式的解，要获得空间中的泄漏气体的浓度需要利用计算机进行数值求解。

步骤S41：将步骤S3建立的非线性偏微分方程（9）转换为用空间离散坐标和时间坐标表示的方程，则可写出向量形式微分方程的差分形式。在将微分方程转换为差分方程时，差分格式的选择非常重要。此处以一阶向前差分为例来描述式（9）的差分方程：

$\frac{U_{i,j,k}^{t+\mathrm{\Δt}}-U_{i,j,k}^t}{\mathrm{\Δt}}=J-\frac{F_{i+1,j,k}^t-F_{i,j,k}^t}{x_{i+1}-x_i}-\frac{G_{i+1,j,k}^t-G_{i,j,k}^t}{y_{i+1}-y_i}-\frac{H_{i+1,j,k}^t-H_{i,j,k}^t}{z_{i+1}-z_i}---\left(10\right).$

步骤S42：给出方程定解条件。方程组的定解条件包括初始条件和边界条件，其中：

若所述泄漏源为连续泄漏源，则初始条件为：t＝0时，泄漏面（在连续泄漏时，把泄漏源视为坐标系中一个泄漏气体的泄漏面，泄漏面的大小和形状可根据需要设置）上泄漏气体扩散浓度为100%，泄漏质量流率为q，方向垂直于泄漏面；固定区计算域内其余各处（除泄漏面之外的区域）空气浓度为100%，速度为V（m/s）（即风速），p取当地大气压值（Pa），ρ为空气密度（Kg/m³）；

在固定区计算域的边界上需要给出边界条件，本实施方式中的边界条件为远场边界条件。远场边界条件可从流体力学的相关资料（例如教材）中查到，此处其值可从大气的气象参数中获得，主要包括风速V（m/s）、当地大气压值p（Pa）和空气密度ρ（Kg/m³）；

若所述泄漏源为瞬时泄漏源，则初始条件为：t＝0时，以泄漏源为中心，体积为的球形区域内，其中ρ对应组分的气体密度，即ρ表示泄漏气体密度，泄漏气体扩散浓度为100%；固定区计算域内其余各处（除该球形区域之外的区域）的值与连续泄漏时相同；

瞬时泄漏时的边界条件设置可以与上述连续泄漏时相同。

步骤S43：根据定解条件迭代求解所述固定区计算域内的泄漏气体扩散浓度分布场的数值解，该数值解中包含参数t（时间坐标）。。

步骤S44：输出步骤S43中得到的泄漏气体扩散浓度分布场的数值解。

步骤S5：建立泄漏气体扩散浓度分布场与采集到的参数中的部分或全部间的函数关系：

若泄漏源为瞬时泄漏源，则建立如式（11）的泄漏气体扩散浓度分布场与x,y,z,t,ρ₀,V,T,H₀,ρ₁,α,Q,L之间的函数关系：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=f(x,y,z,t,{\mathrm{\ρ}}_0,V,T,H_0,{\mathrm{\ρ}}_1,\mathrm{\α},Q,L);---\left(11\right)$

若泄漏源为连续泄漏源，则建立如式（12）的泄漏气体扩散浓度分布场与x,y,z,t,ρ₀,q,V,T,H₀,ρ₁,α,T_t0,L之间的函数关系：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=f(x,y,z,t,{\mathrm{\ρ}}_0,q,V,T,H_0,{\mathrm{\ρ}}_1,\mathrm{\α},T_{t0},L);---\left(12\right)$

其中，ρ₀为泄漏气体密度（Kg/m³），H₀为泄漏源海拔（m），q为泄漏质量流率（Kg/s），Q为泄漏量（Kg/s），T_t0为泄漏时间（s），ρ₁为云层密度（Kg/m³），V为风速（m/s），α为风向（可为风向与正东向之间的夹角，单位为弧度），T为温度（℃），L为泄漏源特征尺寸（m），而f用于表示一种函数关系，并不代表具体的函数形式。

步骤S6：将步骤S5中得到的函数关系表示为无量纲化函数关系，包括：

步骤S61：对式（11）和式（12）进行无量纲化函数关系表示：

若泄漏源为瞬时泄漏源，则建立如式（13）的无量纲化函数关系：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=f(\frac{x}{H_0},\frac{y}{H_0},\frac{z}{H_0},\frac{V}{H_0/t},\frac{{\mathrm{\ρ}}_0}{{\mathrm{QH}}_0^{-3}},\frac{{\mathrm{\ρ}}_1}{{\mathrm{QH}}_0^{-3}},\frac{T}{{\mathrm{QH}}_{0}t},\frac{L}{H_0},\mathrm{\α})---\left(13\right)$

若泄漏源为连续泄漏源，则建立如式（14）的无量纲化函数关系：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=f(\frac{x}{H_0},\frac{y}{H_0},\frac{z}{H_0},\frac{V}{H_0/t},\frac{{\mathrm{\ρ}}_0}{{\mathrm{qtH}}_0^{-3}},\frac{T}{{q{t}^{2}H}_0},\frac{{\mathrm{\ρ}}_1}{\mathrm{qt}{H}_0^{-3}},\frac{T_{t0}}{t},\frac{L}{H_0},\mathrm{\α})---\left(14\right)$

需要说明的是，由于本实施方式中，在步骤S2中设定坐标原点为泄漏源的中心，所以在式（11）～式（14）的方程中未出现泄漏源的经度和纬度。而泄漏源的高度反应出泄漏源与地面间的距离，该距离可影响泄漏气体的扩散过程，所以海拔H₀（单位为m）被体现在方程中。

步骤S62：建立泄漏气体扩散浓度分布场的无量纲化函数关系。

若泄漏源为瞬时泄漏源，则有：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=a_{10}+a_{11}\left(\frac{x}{H_0}\right)+a_{12}{\left(\frac{x}{H_0}\right)}^2+...+a_{{1k}_1}{\left(\frac{x}{H_0}\right)}^{k_1}+...+$

$a_{21}\left(\frac{y}{H_0}\right)+a_{22}{\left(\frac{y}{H_0}\right)}^2+...+a_{{2k}_2}{\left(\frac{y}{H_0}\right)}^{k_2}+...+$

$a_{31}\left(\frac{z}{H_0}\right)+a_{32}{\left(\frac{z}{H_0}\right)}^2+...+a_{{3k}_3}{\left(\frac{z}{H_0}\right)}^{k_3}+...+---\left(15\right)$

$...+$

$a_{91}\left(\mathrm{\α}\right)+a_{92}{\left(\mathrm{\α}\right)}^2+...+a_{{9k}_9}{\left(\mathrm{\α}\right)}^{k_9}$

若泄漏源为连续泄漏源，则有：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=a_{10}+a_{11}\left(\frac{x}{H_0}\right)+a_{12}{\left(\frac{x}{H_0}\right)}^2+...+a_{{1k}_1}{\left(\frac{x}{H_0}\right)}^{k_1}+...+$

$a_{21}\left(\frac{y}{H_0}\right)+a_{22}{\left(\frac{y}{H_0}\right)}^2+...+a_{{2k}_2}{\left(\frac{y}{H_0}\right)}^{k_2}+...+$

$a_{31}\left(\frac{z}{H_0}\right)+a_{32}{\left(\frac{z}{H_0}\right)}^2+...+a_{{3k}_3}{\left(\frac{z}{H_0}\right)}^{k_3}+...+---\left(16\right)$

$...+$

$a_{101}\left(\mathrm{\α}\right)+a_{102}{\left(\mathrm{\α}\right)}^2+...+a_{{10k}_{10}}{\left(\mathrm{\α}\right)}^{k_{10}}$

上式中，在瞬时泄漏的情况下，该无量纲函数关系中共有(k₁+k₂+k₃+...+k₉+1)个待定系数在连续泄漏的情况下，该无量纲函数关系中共有(k₁+k₂+k₃+...+k₁₀+1)个待定系数根据经验，通过权衡结果精度和计算量，可选择k₁、k₂......k₁₀小于等于10。

步骤S63：确定待定系数的值。

当瞬时泄漏时，从步骤S4中计算得到的泄漏气体扩散浓度分布场的数值解中，选取(k₁+k₂+k₃+...+k₉+1)个数值解来建立线性方程组，以求解方程系数

当连续泄漏时，从步骤S4中计算得到的泄漏气体扩散浓度分布场的数值解中，选取(k₁+k₂+k₃+...+k₁₀+1)个数值解来建立线性方程组，以求解方程系数

因此，瞬时泄漏时，有(k₁+k₂+k₃+...+k₉+1)已知值，且对应于每个已知浓度值， $(\frac{x}{H_0},\frac{y}{H_0},\frac{z}{H_0},\frac{V}{H_0/t},\frac{{\mathrm{\ρ}}_0}{{\mathrm{QH}}_0^{-3}},\frac{{\mathrm{\ρ}}_1}{{\mathrm{QH}}_0^{-3}},\frac{T}{{\mathrm{QH}}_{0}t},\frac{L}{H_0},\mathrm{\α})$ 都是唯一的，所以可线性求解系数类似的，对于连续泄漏，有(k₁+k₂+k₃+...+k₁₀+1)已知浓度值，且对应于每个已知值， $(\frac{x}{H_0},\frac{y}{H_0},\frac{z}{H_0},\frac{V}{H_0/t},\frac{{\mathrm{\ρ}}_0}{{\mathrm{qtH}}_0^{-3}},\frac{T}{{q{t}^{2}H}_0},\frac{{\mathrm{\ρ}}_1}{\mathrm{qt}{H}_0^{-3}},\frac{T_{t0}}{t},\frac{L}{H_0},\mathrm{\α})$ 都是唯一的，所以可线性求解系数

步骤S7：将感兴趣的点的空间坐标(x,y,z)和时间t坐标代入步骤S6中得到的所述系数已被确定的所述无量纲化函数关系，计算得到所述感兴趣的点的泄漏气体扩散浓度。因为步骤S6中建立的方程是适用于开放的三维空间的，所以该感兴趣的点可以是固定区计算域内的点或者是固定区计算域外的点。

由此，可以获得空间任意感兴趣的位置的泄漏气体扩散浓度。

步骤S8：计算存在多个泄漏源时的泄漏气体扩散浓度或泄漏气体扩散浓度分布场。

一般来讲，若同时存在多个泄漏源（总数记为B个），且多个泄漏源之间相互独立，可对多个泄漏源中的每个泄漏源分别执行步骤S1至S7，得到感兴趣的点对应于每个泄漏源的泄漏气体扩散浓度，再将所有泄漏源的泄漏气体扩散浓度相加，以得到感兴趣的点的总的泄漏气体扩散浓度；或者针对多个泄漏源中的每一个，按照步骤1～步骤6得到对应于该泄漏源的泄漏气体扩散浓度分布场的无量纲化函数关系，总的泄漏气体扩散浓度分布场是多个泄漏源的泄漏气体扩散浓度分布场的加和，再将所述感兴趣的点的空间坐标和时间坐标代入总的泄漏气体扩散浓度分布场的无量纲化函数关系，得到感兴趣的点的总的泄漏气体扩散浓度。设(X_b,Y_b,H_b)表示第b个泄漏源的经度、纬度和海拔。则总的泄漏气体扩散浓度场 $E\left\{M_{i_{\mathrm{Atm}}}(x,y,z,t)\right\}$ 为：

$E\left\{M_{i_{\mathrm{Atm}}}(x,y,z,t)\right\}=\underset{b=0,...,B-1}{\mathrm{\Σ}}{M}_{i_{\mathrm{Atm}}}{(x,y,z,t)}_{(X_b,Y_b,H_b)}---\left(17\right)$

要注意的是，在对多个泄漏源的泄漏气体扩散浓度分布场进行加和时，如果针对不同泄漏源选择的坐标系不同，例如所选择的坐标原点不同，应该先进行坐标系转换，得到同一坐标系下各泄漏源的泄漏气体扩散浓度场的无量纲化函数关系，再将各个浓度场相加。

本发明所述的方法与现有技术相比，在针对单一泄漏源进行计算时得到的结果更加可靠和健壮，因此，对多泄漏源进行叠加求解时与现有技术相比能得到的更为准确的结果。本发明所述的方法已经在上述的计算机及软件平台上进行了方法的有效性得到了验证,并以地理信息系统（GIS）获得参数作为输入进行了方法验证，目前已部署在相关用户单位。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种用于三维空间的泄漏气体扩散预测方法，该方法包括：

步骤S1，采集参数，所述参数包括地理环境参数和泄漏源自身的参数；

步骤S2，建立坐标系，并参照该坐标系建立固定区计算域，将所述固定区计算域内的气体离散化为由多个离散点界定的流体微团；

步骤S3，基于流体力学领域的纳维-斯托克斯方程建立在所述固定区计算域内的泄漏气体扩散方程，所述泄漏气体扩散方程中包括质量守恒方程、动量方程和能量守恒方程；

步骤S4，将所述步骤S3中得到的所述质量守恒方程、所述动量方程和所述能量守恒方程转换为关于空间离散坐标和时间坐标的差分方程，并根据定解条件和采集到的所述参数，以迭代的方式求解所述固定区计算域内的泄漏气体扩散浓度分布场的数值解；

步骤S5，建立所述泄漏气体扩散浓度分布场与采集到的所述参数中的部分或全部间的函数关系；

步骤S6，将所述步骤S5中得到的所述函数关系表示为包含待定系数的无量纲化函数关系，并且将从所述步骤S4中得到的所述数值解中的部分或全部代入所述无量纲化函数关系，以确定所述待定系数的值；以及

步骤S7，将感兴趣的点的空间坐标和时间坐标代入步骤S6中得到的所述系数已被确定的所述无量纲化函数关系，计算得到所述感兴趣的点的泄漏气体扩散浓度，所述感兴趣的点在所述固定区计算域内或者在所述固定区计算域外。

2.根据权利要求1所述的泄漏气体扩散预测方法，还包括：

当所述泄漏源为多个泄漏源时，针对所述多个泄漏源中的每一个，按照所述步骤1～步骤7计算该泄漏源在所述感兴趣的点的所述泄漏气体扩散浓度，将计算得到的所述多个泄漏源的所述泄漏气体扩散浓度相加，以得到所述感兴趣的点的总的泄漏气体扩散浓度；或者针对所述多个泄漏源中的每一个，按照所述步骤1～步骤6得到对应于该泄漏源的所述泄漏气体扩散浓度分布场的所述无量纲化函数关系，总的泄漏气体扩散浓度分布场是所述多个泄漏源的所述泄漏气体扩散浓度分布场的和，再将所述感兴趣的点的空间坐标和时间坐标代入所述总的泄漏气体扩散浓度分布场的所述无量纲化函数关系，以得到所述感兴趣的点的总的泄漏气体扩散浓度。

3.根据权利要求1所述的泄漏气体扩散预测方法，其中所述地理环境参数包括泄漏源所在地的经度X₀、纬度Y₀和海拔H₀、云层密度ρ₁、风速V、风向α和温度T；所述泄漏源自身的参数包括泄漏气体密度ρ₀和泄漏源特征尺寸L，如果所述泄漏源是瞬时泄漏源，还包括泄漏量Q，如果所述泄漏源是连续泄漏源，还包括泄漏质量流率q和泄漏时间T_t0。

4.根据权利要求3所述的泄漏气体扩散预测方法，其中在所述步骤3中建立所述固定区计算域内的所述泄漏气体扩散方程时，对所述固定区计算域内的所述泄漏气体做层流假设，并假设所述流体微团在运动过程中是绝能的。

5.根据权利要求3所述的泄漏气体扩散预测方法，其中所述步骤2包括：

步骤S21：以所述泄漏源的中心为坐标原点O，以所述风向作为直角坐标系的X正方向，建立三维右手坐标系OXYZ；

步骤S22：以所述泄漏源特征尺寸L为基准，建立以所述坐标原点O为中心的正方体计算域作为所述固定区计算域，坐标原点O到正方体各平面之间的垂直距离为βL，其中，β为远大于1，且小于等于1000的一个数值，相应的计算域为(2βL)×(2βL)×(2βL)的正方体；

步骤S23：对所述步骤S22得到的所述正方体计算域沿X轴、Y轴和Z轴进行离散化，离散化后的间隔分别为Δx，Δy和Δz，相应的X轴、Y轴和Z轴分别对应的点数为 和所述离散点的计算域空间坐标为(x,y,z)，其中 $x=\mathrm{\Δx}\×(-\frac{N_x}{2}:1:(\frac{N_x}{2}-1)),$ $y=\mathrm{\Δy}\×(-\frac{N_y}{2}:1:(\frac{N_y}{2}-1)),z=\mathrm{\Δz}\×(-\frac{N_z}{2}:1:(\frac{N_z}{2}-1)).$

6.根据权利要求4所述的泄漏气体扩散预测方法，其中所述步骤S3建立的所述泄漏气体方程包括：

$J_{i_{\mathrm{Atm}}}=-D_{i_{\mathrm{Atm}},m}\▿M_{i_{\mathrm{Atm}}};$

方程7：p＝ρRT;以及

e＝c_vT;

方程8： $\frac{\∂U}{\∂t}=J-\frac{\∂F}{\∂x}-\frac{\∂G}{\∂y}-\frac{\∂H}{\∂z}$

其中， $U=\left\{\begin{array}{c}{\mathrm{\ρM}}_{i_{\mathrm{Atm}}}\\ \mathrm{\ρu}\\ \mathrm{\ρv}\\ \mathrm{\ρw}\\ \mathrm{\ρ}(e+\frac{V^2}{2})\end{array}\right\},$ $F=\left\{\begin{array}{c}{\mathrm{\ρM}}_{i_{\mathrm{Atm}}}u+J_{i_{\mathrm{Am}}}^x\\ {\mathrm{\ρu}}^2+p\\ \mathrm{\ρuv}\\ \mathrm{\ρuw}\\ \mathrm{\ρu}(e+\frac{V^2}{2})+\mathrm{pu}\end{array}\right\},$ $F=\left\{\begin{array}{c}{\mathrm{\ρM}}_{i_{\mathrm{Atm}}}v+J_{i_{\mathrm{Am}}}^y\\ \mathrm{\ρvu}\\ {\mathrm{\ρv}}^2+p\\ \mathrm{\ρvw}\\ \mathrm{\ρv}(e+\frac{V^2}{2})+\mathrm{pv}\end{array}\right\},$ $H=\left\{\begin{array}{c}{\mathrm{\ρM}}_{i_{\mathrm{Atm}}}w+J_{i_{\mathrm{Am}}}^z\\ \mathrm{\ρwu}\\ \mathrm{\ρwv}\\ {\mathrm{\ρw}}^2+p\\ \mathrm{\ρw}(e+\frac{V^2}{2})+\mathrm{pw}\end{array}\right\},$ $J=\left\{\begin{array}{c}0\\ {\mathrm{\ρf}}_x\\ {\mathrm{\ρf}}_y\\ {\mathrm{\ρf}}_z\\ \mathrm{\ρ}({\mathrm{uf}}_x+{\mathrm{vf}}_y+{\mathrm{wf}}_z)\end{array}\right\}.$

7.根据权利要求1所述的泄漏气体扩散预测方法，其中所述步骤S4中所述差分方程采用一阶向前差分的差分格式。

8.根据权利要求6所述的泄漏气体扩散预测方法，其中所述步骤S4包括：

步骤S41：将所述步骤S3中得到的所述方程8转换为关于空间离散坐标和时间坐标的一阶向前差分方程：

$\frac{U_{i,j,k}^{t+\mathrm{\Δt}}-U_{i,j,k}^t}{\mathrm{\Δt}}=J-\frac{F_{i+1,j,k}^t-F_{i,j,k}^t}{x_{i+1}-x_i}-\frac{G_{i+1,j,k}^t-G_{i,j,k}^t}{y_{i+1}-y_i}-\frac{H_{i+1,j,k}^t-H_{i,j,k}^t}{z_{i+1}-z_i};$

步骤S42：给出方程定解条件，所述定解条件包括初始条件和边界条件，所述边界条件为远场边界条件，

若所述泄漏源为连续泄漏源，则所述初始条件包括：t＝0，泄漏面上的泄漏气体扩散浓度为100%，泄漏质量流率为采集到的所述泄漏质量流率q，泄漏方向垂直于所述泄漏面；所述固定区计算域内其余各处的空气浓度为100%，所述其余各处的气体运动速度为所述风速V，所述其余各处的压强为当地大气压值p，所述其余各处的气体密度为空气密度ρ；

若所述泄漏源为瞬时泄漏源，则所述初始条件包括，t＝0，以所述坐标原点为中心，体积为所述泄漏量Q除以所述泄漏气体密度ρ₀的球形区域内，泄漏气体扩散浓度为100%；所述固定区计算域内的所述其余各处的所述初始条件和连续泄漏时相同；

步骤S43：根据所述定解条件迭代求解所述固定区计算域内的所述泄漏气体扩散浓度分布场的所述数值解，该数值解中包含时间坐标t。

9.根据权利要求3所述的泄漏气体扩散预测方法，其中，所述步骤S5中建立的所述函数关系为：

若所述泄漏源为瞬时泄漏源，则所述函数关系为：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=f(x,y,z,t,{\mathrm{\ρ}}_0,V,T,H_0,{\mathrm{\ρ}}_1,\mathrm{\α},Q,L);$

若所述泄漏源为连续泄漏源，则所述函数关系为：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=f(x,y,z,t,{\mathrm{\ρ}}_0,q,V,T,H_0,{\mathrm{\ρ}}_1,\mathrm{\α},T_{t0},L);$ 以及

其中，是所述泄漏气体扩散浓度分布场。

10.根据权利要求9所述的泄漏气体扩散预测方法，其中所述步骤S6包括：

步骤S61：将所述步骤S5中得到的所述函数关系进行无量纲化表示：

若所述泄漏源为瞬时泄漏源，则所述无量纲化函数关系为：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=f(\frac{x}{H_0},\frac{y}{H_0},\frac{z}{H_0},\frac{V}{H_0/t},\frac{{\mathrm{\ρ}}_0}{{\mathrm{QH}}_0^{-3}},\frac{{\mathrm{\ρ}}_1}{{\mathrm{QH}}_0^{-3}},\frac{T}{{\mathrm{QH}}_{0}t},\frac{L}{H_0},\mathrm{\α}),$

若所述泄漏源为连续泄漏源，则所述无量纲化函数关系为：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=f(\frac{x}{H_0},\frac{y}{H_0},\frac{z}{H_0},\frac{V}{H_0/t},\frac{{\mathrm{\ρ}}_0}{{\mathrm{qtH}}_0^{-3}},\frac{T}{{q{t}^{2}H}_0},\frac{{\mathrm{\ρ}}_1}{\mathrm{qt}{H}_0^{-3}},\frac{T_{t0}}{t},\frac{L}{H_0},\mathrm{\α});$

步骤S62：建立泄漏气体扩散浓度分布场的无量纲化函数关系：

若所述泄漏源为瞬时泄漏源，则：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=a_{10}+a_{11}\left(\frac{x}{H_0}\right)+a_{12}{\left(\frac{x}{H_0}\right)}^2+...+a_{{1k}_1}{\left(\frac{x}{H_0}\right)}^{k_1}+...+$

$a_{21}\left(\frac{y}{H_0}\right)+a_{22}{\left(\frac{y}{H_0}\right)}^2+...+a_{{2k}_2}{\left(\frac{y}{H_0}\right)}^{k_2}+...+$

$a_{31}\left(\frac{z}{H_0}\right)+a_{32}{\left(\frac{z}{H_0}\right)}^2+...+a_{{3k}_3}{\left(\frac{z}{H_0}\right)}^{k_3}+...+,$

$...+$

$a_{91}\left(\mathrm{\α}\right)+a_{92}{\left(\mathrm{\α}\right)}^2+...+a_{{9k}_9}{\left(\mathrm{\α}\right)}^{k_9}$

若所述泄漏源为连续泄漏源，则：

$M_{i_{\mathrm{Atm}}}(x,y,z,t)=a_{10}+a_{11}\left(\frac{x}{H_0}\right)+a_{12}{\left(\frac{x}{H_0}\right)}^2+...+a_{{1k}_1}{\left(\frac{x}{H_0}\right)}^{k_1}+...+$

$a_{21}\left(\frac{y}{H_0}\right)+a_{22}{\left(\frac{y}{H_0}\right)}^2+...+a_{{2k}_2}{\left(\frac{y}{H_0}\right)}^{k_2}+...+$

$a_{31}\left(\frac{z}{H_0}\right)+a_{32}{\left(\frac{z}{H_0}\right)}^2+...+a_{{3k}_3}{\left(\frac{z}{H_0}\right)}^{k_3}+...+.$

$...+$

$a_{101}\left(\mathrm{\α}\right)+a_{102}{\left(\mathrm{\α}\right)}^2+...+a_{{10k}_{10}}{\left(\mathrm{\α}\right)}^{k_{10}}$

其中，或是待定系数；

步骤S63：确定所述待定系数；

若所述泄漏源为瞬时泄漏源，从所述步骤S4中得到的所述数值解中选取(k₁+k₂+k₃+...+k₉+1)个所述数值解，代入步骤S62中得到所述无量纲化函数关系，求解

若所述泄漏源为连续泄漏源，从所述步骤S4中得到的所述数值解中选取(k₁+k₂+k₃+...+k₁₀+1)个所述所述数值解，代入步骤S62中得到的所述无量纲化函数关系，求解