CN106650120A - 一种液烃管道泄漏量测算方法 - Google Patents

Abstract

本发明涉及一种液烃管道泄漏量测算方法，特征在于其包括以下步骤：1)建立泄漏管道的通用物理模型；2)建立泄漏管道中各管段的连续数学模型，得到描述管内油流的控制方程；3)对泄漏管道和时间域进行网格划分得到离散的计算区域；4)将连续数学模型进行离散化，得到描述管内油流控制方程的特征方程；5)根据网格划分结果，将管内油流控制方程的特征方程转换为离散格式；6)根据通用物理模型中各管段的相邻处不同的边界类型，确定不同边界条件并建立其离散格式；7)根据管内油流控制方程的离散格式以及各边界条件对管道全线进行水热力耦合模拟计算，得到泄漏管道的泄漏量。本发明可以广泛应用于液烃管道泄漏量的测算领域。

Description

一种液烃管道泄漏量测算方法

技术领域

本发明涉及管道泄漏领域，特别是关于一种液烃管道泄漏量测算方法。

背景技术

能源行业输送液态流体的管道主要包括原油管道、成品油管道、集输管道、轻烃管道、LPG管道和LNG管道等，均属于液态烃类管道。液烃管道因其地位特殊，管中液体极易燃易爆，又常大流量高压输送，泄漏后果严重，潜在危害大。国内外与管道泄漏相关的研究主要集中在泄漏检测与定位领域，针对泄漏量计算的相关研究几乎是空白的。同时，针对管道泄漏过程的相关研究也只停留在较为基础的阶段，因考虑的影响因素较少、对泄漏过程流体流动状态刻画较为简单，研究结果与实际管道的泄漏过程相差较大。

目前研究液烃在泄漏过程中的累计泄漏量方法较少，仅有的相关研究提出可通过稳态计算、数值模拟、统计分析、事后估计、实验分析等泄漏量计算方法，但都存在误差较大或无法应用于工程实际的缺点。下面对几种主要的计算方法做简单介绍。

(1)解析法

将泄漏过程当作一维稳定管流的小孔出流，利用伯努利方程计算水平和倾斜管道液烃在泄漏点的泄漏速率与泄漏量。解析法求解泄漏量时只能计算出管道处于稳定运行状态下的泄漏总量，计算方法简单，但是误差大，且该方法只有事先确定发生泄漏的时间才能准确计算出泄漏量。

(2)数值模拟法

数值模拟法是基于流体力学的基本理论，建立微分方程模型，利用数值模拟方法计算泄漏量。该方法的不足之处为模型中的环境压力及泄漏孔出流系数常为定值，无法反映孔口出流与管外渗流扩散的耦合效应，也未考虑液烃温度变化对泄漏产生的影响，且模型规模变大之后，数值模拟难度大，如长管道、多边界等很难快速模拟。

(3)统计学法

影响泄漏量的因素多，因此有学者采用信息系统或统计学方法预测泄漏量，如KimB利用统计模型预测管道发生故障时的漏油量。基于特定输油管道的基础参数和故障数据，采用双参数韦伯分布(Wei Bull Distribution)预测液烃的泄漏量。但是该方法必须有足够的样本数据来构造韦伯分布，且只能针对特定管道对泄漏事故中的总泄漏量给出估计，通用性较差。

(4)间接计算法

直接计算难度大，因此有学者提出了间接计算泄漏量的模型，刘国华提出在确定了浸油土壤边际曲线y＝ax²后，将其绕y轴旋转360°后，再由定积分旋转体体积公式求得浸油土壤体积。根据测算的含油土壤体积，通过环形分层多点取样采样方法得到单位浸油土壤的密度及含油率，估算液烃总泄漏量。但是此方法是事后估计法，无法给出泄漏过程中的泄漏速率和实时累计泄漏量，同时，模型也没有考虑泄漏过程中液烃蒸气在空气中的扩散损失，因此实际应用价值不足。

(5)实验分析法

张蓉爱在实验室建立管道泄漏实验模型，确立了泄漏量与泄漏孔的大小、流体的密度、黏度、管内压力以及流体的平均流速等因素的关系，但数学模型是基于伯努利方程的小孔稳态出流模型，实验测定的是出流系数。Jasper A用实验对上、下游阀门操作时的泄漏量的影响进行了研究，泄漏量靠重力仪测得，并未给出理论分析和计算模型。付建民搭建液相管道小孔泄漏实验系统研究管道流量及压力对管内液体压力及泄漏速率变化的影响规律，提出了小孔泄漏稳定压力计算经验公式，但并未给出泄漏量计算公式，也不能用于预测泄漏过程瞬态压力变化。然而，实验研究无法真实还原管道泄漏的实际工况，且存在多种人为误差，无法用于实际管道泄漏量的测算。

发明内容

针对上述问题，本发明的目的是提供一种液烃管道泄漏量测算方法，针对液烃管道泄漏过程中的泄漏量进行实时测算，分析了多种边界条件下的边界点计算方法，准确性高。

为实现上述目的，本发明采取以下技术方案：一种液烃管道泄漏量测算方法，特征在于其包括以下步骤：1)首先建立泄漏管道的通用物理模型；2)根据建立的泄漏管道的通用物理模型，建立泄漏管道中各管段的连续数学模型，得到描述管内油流的控制方程；3)采用各管段统一时步法对泄漏管道和时间域进行网格划分，得到离散的计算区域；4)采用特征线法，将步骤2)中建立的泄漏管道的连续数学模型进行离散化，得到描述管内油流控制方程的特征方程；5)根据步骤3)中的网格划分结果，将步骤4)中得到的管内油流控制方程的特征方程转换为离散格式；6)根据步骤1)建立的通用物理模型中各管段的相邻处不同的边界类型，确定不同边界条件，并建立其离散格式；7)根据步骤5)中得到的管内油流控制方程的离散格式以及步骤6)中确定的各边界条件对管道全线进行水热力耦合模拟计算，得到泄漏过程中所需相应计算时间内泄漏管道的泄漏量。

所述步骤2)中，管内油流的控制方程包括管内油流的连续性方程、动量方程和换热方程；其中，管内油流的连续性方程为：

管内油流的动量方程为：

管内油流的换热方程为：

式中，t为时间，单位为s；ρ为油品在管截面上的平均密度，单位为kg/m³；x为距泄漏管段起点的距离，单位为m；v为管内油品的平均速度，单位为m/s；A为泄漏管段截面积，单位为m²；g为重力加速度，单位为m/s²；θ是泄漏管段与水平方向的夹角，单位为rad；D为管内径，单位为m；λ为达西摩阻系数；q为油品与单位面积管壁单位时间内的热流密度，单位为W/m²，T为管内油品平均温度，单位为℃；α_p为油品的膨胀系数，单位为℃^-1；c为油品的比热容，单位为J/(kg·℃)；a_o,g为压力波在不同介质中的传播速度，单位为m/s；ρ_o为不同批次的油品密度，单位为kg/m³；k_o为对应批次油品的弹性模量，单位为Pa；D_g为不同管段直径，单位为m；E_g为对应管段的杨氏弹性模量，单位为Pa；δ_g为对应管段的壁厚，单位为m。

所述步骤4)中，得到的管内油流控制方程的特征方程为：

C⁺特征方程：

C^-特征方程：

v特征方程：

所述步骤5)中，得到的所述控制方程的离散格式为：

C⁺特征方程：

C^-特征方程：

v特征方程：

式中，Δx为空间步长，单位为m；Δt为时间步长，单位为s；f为列宾宗摩阻系数；w分别为管道横截面，单位为m²；分别为节点P处j+1时刻的流量，压头，温度和压力；分别为节点A处j时刻的流量，压头，温度和压力；分别为节点B处j时刻的流量和压头。

所述步骤6)中，得到的边界条件包括上下游边界、阀门边界、变径点边界、泄漏点边界、混油边界和液柱分离边界。

所述上下游边界中，上、下游边界节点的压头和流量计算公式分别为：

式中，分别为管段i的第N个节点处j+1时刻的压头和流量，也即管段i上游边界的压头和流量；分别为管段i的第N-1个节点处j时刻的压头和流量方程的系数；分别为管段i+1的第0个节点处j+1时刻的压头和流量，也即管段i上游边界的压头和流量；分别为管段i+1的第1个节点处j时刻的压头和流量方程的系数；

所述阀门边界中，阀门边界节点的压头和流量的计算公式分别为：

式中，分别为阀门上游边界在j+1时刻的压头和流量；分别为阀门上游边界j时刻的压头和流量方程的系数；分别为阀门下游边界j+1时刻的压头和流量；分别为阀门下游边界j时刻的压头和流量方程的系数；Q_p为阀门处的流量；

所述变径点边界中，变径点处压头和流量的计算公式为：

式中，分别为变径点上游边界j+1时刻的压头和流量；分别为阀门上游边界j时刻的压头和流量方程的系数；分别为阀门下游边界j+1时刻的压头和流量；分别为阀门下游边界j时刻的压头和流量方程的系数；H_p、Q_p分别为变径点上下游处的压头和流量；

所述泄漏点边界中，泄漏点处压头和流量的计算公式为：

且泄漏点处压头H_p为：

式中， 分别为泄漏点上游边界j+1时刻的压头和流量；分别为泄漏点上游边界j时刻的压头和流量方程的系数；分别为泄漏点下游边界j+1时刻的压头和流量；分别为泄漏点下游边界j时刻的压头和流量方程的系数；H₀为泄漏孔外部的压头，单位为m；C₀为泄漏孔的流量系数，且C₀＝αC_dw≈0.6～0.65w；q_P为泄漏孔处瞬时泄漏流量，单位为m³/s；

所述混油边界中，混油边界的压力和流量计算公式为：

其中，ρ_A为前行油品的密度，kg/m³；ρ_B为后行油品的密度，kg/m³；β＝ρ_A/ρ_B；Z为混油边界处的高程，单位为m；分别为混油界面点上游边界j+1时刻的压头和流量；分别为混油界面点上游边界j时刻的压头和流量方程的系数；分别为混油界面点下游边界j+1时刻的压头和流量；分别为混油界面点下游边界j时刻的压头和流量方程的系数；P_p、Q_p分别为混油界面点的压力和流量；

所述液柱分离边界中，液柱分离边界节点的计算方法为：

式中，分别为液柱分离点上游边界在j+1时刻的压头和流量；分别为泄漏点下游边界在j+1时刻的压头和流量；分别为液柱分离点上、下游边界在j时刻的流量；H_vapo_r为液体饱和蒸汽压力，单位为m；分别为液柱分离点在j时刻的速度；ψ为加权因子，且0.5<ψ<1。

所述步骤7)中，对管道全线进行水热力耦合模拟计算的方法包括以下步骤：①首先确定管道全线的基本参数，包括管道直径、管长、管道里程和高程；②根据管道全线的基本参数，确定该管道物理模型中的相应物理单元，根据步骤6)中的边界条件处理方法，设置相应的边界条件和初始条件；③根据当前温度，计算管道内油品的物性参数，包括油品密度ρ、粘度μ和列宾宗摩阻系数f；④判断当前节点是否发生液柱分离，若发生液柱分离，则采用液柱分离的边界条件计算该节点当前时步的压力和流量值；若没有发生液柱分离，则根据当前时步各节点的油品物性，计算管道各节点当前时步的压力和流量值；⑤根据各节点当前时步的压力和流量值，采用v特征方程计算各节点油品当前时步的温度值；⑥判断所得到的温度值是否收敛，如果是，表明当前时步下各管段节点的温度均计算完毕，进入步骤⑦，如果否，则返回步骤③，将得到的温度值返回继续进行油品物性计算，直到收敛；⑦确定当前时步油品批次界面位置，并根据得到的温度值以及泄漏点的边界条件计算当前时步的泄漏量；⑧对以后每一时步重复步骤③～⑦，直至达到设定的计算时间。

本发明由于采取以上技术方案，其具有以下优点：1、本发明由于建立了泄漏管道的通用物理模型，包括了各种管道元件对泄漏管道的影响，可以准确描述管道破损处流体的流动状态。2、本发明由于根据物理模型对泄漏出流瞬态过程进行研究，建立了泄漏过程的数学模型，得到描述管内油流的控制方程，实现了对泄漏过程中的泄漏状态进行描述。3、本发明由于针对泄漏管道中不同管段的边界条件进行了研究，给出了不同边界条件下的计算方法，使得泄漏量的计算更加精确。4、本发明由于将泄漏管道处的控制方程和各种边界条件结合起来，对管道全线进行水热力耦合模拟，得到了不同状态下管道泄漏量的实时计算，计算结果更加真实可靠。本发明可以广泛应用于液烃管道泄漏测算领域中。

附图说明

图1是本发明通用物理模型；

图2是本发明特征线与有限差分法；

图3是本发明内部节点示意图；

图4是本发明上游边界节点计算方法示意图；

图5是本发明下游边界节点计算方法示意图；

图6是本发明阀门边界节点计算方法示意图；

图7是本发明变径点边界节点计算方法示意图；

图8是本发明泄漏边界节点计算方法示意图；

图9是本发明液柱分离边界节点计算方法示意图；

图10是本发明混油边界节点计算方法示意图；

图11是本发明管道泄漏量测算方法流程图；

图12是本发明实施例中贵阳-安顺管道纵断面图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明提出一种液烃管道泄漏量测算方法，包括以下步骤：

1)首先建立泄漏管道的通用物理模型。

管道泄漏发生后管内油品即以泄漏点为源向两端发射压力波。泄漏后对全线管道运行参数的影响是，上下游压力均降低，上游流量增大，下游流量减小。为保证能够得到现场泄漏过程中泄漏管段两端压力(作为边界条件)，选择上、下游泵站间的管道作为计算的基本单元。

如图1所示，为保证物理模型的通用性，泄漏管道中包含的物理单元依次为上游泵站出口1、变径点2、两阀门3、混油界面4、泄漏点5以及下游泵站入口6。各物理单元构成管道内部的多个边界，将泄漏管道相隔为多个管段。

2)根据建立的泄漏管道的通用物理模型，建立泄漏管道中各管段的连续数学模型，得到描述管内油流的控制方程。

泄漏管道的数学模型涉及管内油流的连续性方程、动量方程和能量方程。其中，管内油流的连续性方程为：

管内油流的动量方程为：

管内油流的能量方程为：

q＝-K(T-T₀)(4)

将管内油流的连续性方程、能量方程和动量方程联立，即可得到管内油流的换热方程，管内油流的连续性方程、动量方程和换热方程共同构成管内油流的控制方程。管内油流的换热方程为：

式(1)～式(6)中，t为时间，单位为s；ρ为油品在管截面上的平均密度，单位为kg/m³；x为距泄漏管段起点的距离，单位为m；v为管内油品的平均速度，单位为m/s；A为泄漏管段截面积，单位为m²。g为重力加速度，单位为m/s²；θ是泄漏管段与水平方向的夹角，单位为rad；D为管内径，单位为m；λ为达西摩阻系数；e为油品比内能，单位为J/kg；h为油品比焓，单位为J/kg；s为相邻计算节点间的高程差，单位为m；q为油品与单位面积管壁单位时间内的热流密度，单位为W/m²，T为管内油品平均温度，单位为℃；T₀为外界环境温度，单位为℃；K为总传热系数，单位为W/(m²·℃)。α_p为油品的膨胀系数，单位为℃^-1；c为油品的比热容，单位为J/(kg·℃)；a_o,g为压力波在不同介质中的传播速度，单位为m/s；ρ_o为不同批次的油品密度，单位为kg/m³；k_o为对应批次油品的弹性模量，单位为Pa；D_g为不同管段直径，单位为m；E_g为对应管段的杨氏弹性模量，单位为Pa；δ_g为对应管段的壁厚，单位为m。

3)采用各管段统一时步法对泄漏管道和时间域进行网格划分，得到离散的计算区域。

对于复杂管道，一般采用各管段统一时步法进行计算。由步骤1)建立的物理模型可知，泄漏管道存在多个内部边界，各内部边界把管道相隔为多个管段，每个管段因油品物性、管道参数、油品状态参数不同，各管段的波速不同。设相邻两条管段的长度、分段数、波速分别为l₁、n₁、a₁，l₂、n₂、a₂，根据矩形网格计算的要求，计算距步与波速和选取时间步长有关系，即：

其中，n₁和n₂必须为整数，若采用各管段统一时步法进行计算，还必须满足：

Δt₁＝Δt₂(9)

其中，Δx₁、Δx₂分别是两相邻管段的距步，Δt₁、Δt₂分别是两相邻管段的时间步长。

要同时满足上述式(6)、(7)、(8)三个条件，尤其在管段较多、输送多种油品时，往往需采用修改波速或管长、间距内插、向后时间内插、向前时间内插、两步法等方法进行处理，其中间距内插法、向后时间差分法较为常用。

如图2所示，为得到稳定、准确的数值解，本发明采用间距内插法对时间域及管道全线进行网格划分。其中，横坐标表示管道，共划分为n段，如图中的1、2、i、...、n-1、n，空间步长为Δx。纵坐标表示时间，时间步长为Δt。

4)采用特征线法，将步骤2)中建立的泄漏管道的连续数学模型进行离散，得到描述管内油流控制方程的特征方程。

特征线法(method of characteri st ics)是一种基于特征理论的求解双曲型偏微分方程组的近似计算方法，其具有理论严密，物理意义明确，适用范围广等特点。特征线法针对长输管道系统的各边界节点及内部网格点独立建立相应的低维线性或非线性代数方程组，编程求解较为简单。

而管道内非稳定流动的基本方程为一组拟线性双曲型偏微分方程，非稳定泄漏过程涉及的问题为快瞬变流动问题。而快瞬变流动的特点是在短时间内完成流动状态的转变，且流动参数发生显著变化。考虑到瞬变过程中的水力参数随时间变化很快，故计算过程中应选取较短的时步，便于研究每一时刻流动参数的变化。显式差分方法是求解上述问题的较好算法，同时可保证很好的稳定性和较高的准确度。

因而，本发明采用特征线法与有限差分相结合的方法，对水热力耦合瞬变流动问题进行求解，并对算法的稳定性和计算精度进行检验。根据特征方程一般形式，采用特征线法，将管内油流的连续性方程、动量方程和换热方程转换为描述管内油流的连续性方程、动量方程和换热方程的特征方程。

C⁺特征方程：

C^-特征方程：

v特征方程：

5)根据步骤3)中的网格划分结果，将步骤4)中得到的管内油流控制方程的特征方程转换为离散格式。

如图3所示，为管段内部节点示意图，管道内部介质流动方向如图中黑色箭头线所示。由于各管段内部节点的特点是压力和流量均相同，其中j为时间点，A、B、P为管道内部节点，采用工程上常用的压头H(x,t)、流量Q(x,t)代替压力P(x,t)、速度v(x,t)，从而将特征方程转换为离散格式。其中，P＝ρgH，v＝Q/A。

管内油流的连续性方程、动量方程和换热方程的的离散格式为：

C⁺特征方程：

C^-特征方程：

v特征方程：

式中，Δx为空间步长，单位为m；Δt为时间步长,单位为s；f为列宾宗摩阻系数、w分别为管道横截面，单位为m²；分别为节点P处j+1时刻的流量，压头，温度和压力；分别为节点A处j时刻的流量，压头，温度和压力；分别为节点B处j时刻的流量和压头。

6)根据步骤1)建立的通用物理模型中各管段的相邻处不同的边界类型，确定不同边界条件，并建立其离散格式。

成品油顺序输送水热力耦合数值模拟需基于特定的边界条件。不同的边界类型有其不同的特点。本发明根据步骤1)中建立的物理模型，针对不同的边界类型确定不同的边界条件，并建立其离散格式。主要包括以下几种边界条件：

①上下游边界

如图4、图5所示，为依据通用物理模型建立的上下游边界示意图，其中j为时间点，i为管段，0、1、N-1、N均为节点，得到上、下游边界节点的压头和流量计算公式分别为：

式(16)～(17)中，分别为管段i的第N个节点处j+1时刻的压头和流量，也即管段i上游边界的压头和流量；分别为管段i的第N-1个节点处j时刻的压头和流量方程的系数；分别为管段i+1的第0个节点处j+1时刻的压头和流量，也即管段i上游边界的压头和流量；分别为管段i+1的第1个节点处j时刻的压头和流量方程的系数。

②阀门边界

如图6所示，为阀门边界示意图。阀门边界又称为扰动边界，它本身的特性是随时间改变的，同时，阀门边界的上下游流量值相等，得到阀门边界节点处压头和流量的计算公式分别为：

式(18)～(21)中，分别为阀门上游边界在j+1时刻的压头和流量；分别为阀门上游边界j时刻的压头和流量方程的系数；分别为阀门下游边界在j+1时刻的压头和流量；分别为阀门下游边界j时刻的压头和流量方程的系数；Q_p为阀门处的流量。

③变径点边界

如图7所示，因变径点两侧管径不同，该点两侧的水力特征线不同，压力波在变径点处发生反射。变径点边界的特点是，边界上下游的压力和流量均相同，因而变径点处压头和流量的计算公式为：

式(22)～(25)中，分别为变径点上游边界j+1时刻的压头和流量；分别为阀门上游边界j时刻的压头和流量方程的系数；分别为阀门下游边界j+1时刻的压头和流量；分别为阀门下游边界j时刻的压头和流量方程的系数；H_p、Q_p分别为变径点上下游处的压头和流量。

④泄漏点边界

假定泄漏具有孔口出流的性质，则有：

式(26)中，H₀为泄漏孔外部的压头，单位为m；C₀为泄漏孔的流量系数，且C₀＝αC_dw≈0.6～0.65w；α为流束收缩系数，其值为0.62～0.66；C_d为孔口流速系数，值为0.98～0.99；w为泄漏孔面积，单位为m²，q_P为泄漏孔处瞬时泄漏流量，m³/s。

如图8所示，为泄漏点边界示意图，泄漏点处压头和流量的计算公式为：

根据式(29)和式(30)可得到泄漏点处压头为：

式中， 分别为泄漏点上游边界j+1时刻的压头和流量；分别为泄漏点上游边界j时刻的压头和流量方程的系数；分别为泄漏点下游边界j+1时刻的压头和流量；分别为泄漏点下游边界j时刻的压头和流量方程的系数。

⑤混油边界

把混油段中央放置于相应的差分节点上，以它作为内部边界，认为节点上游管段内全是后行油品，下游管段全是前行油品。当混油界面移动到下一节点时，混油界面边界相应的移至该节点。因节点之间的长度(即距步)远大于每一时步混油界面移动的距离，故混油界面由当前节点移动至下一节点需要经历很长一段时间，在此不考虑混油界面移动至节点之间时对计算带来的影响。

如图9所示，为混油边界示意图，混油边界的特点是该节点处压力和流量均相同，其压力和流量计算公式为：

其中，ρ_A为前行油品的密度，kg/m³；ρ_B为后行油品的密度，kg/m³；β＝ρ_A/ρ_B；Z为混油边界处的高程，单位为m。分别为混油界面点上游边界j+1时刻的压力和流量；分别为混油界面点上游边界j时刻的压头和流量方程的系数；分别为混油界面点下游边界j+1时刻的压力和流量；分别为混油界面点下游边界j时刻的压头和流量方程的系数；P_p、Q_p分别为混油界面点的压力和流量。

⑥液柱分离边界

除了上述各物理单元构成的各内部边界外，还需要对全线每一点进行液柱分离判别。汽穴流和液柱分离甚为复杂，本发明采用常用的解决方法对液柱分离现象发生时的边界条件进行处理，即只要液体的绝对压力降至其饱和蒸汽压力即发生液柱分离，就把它作为内部边界进行处理。该方法不计及气体的析出，把液柱分离现象中的气体按纯蒸汽处理。

当节点P的绝对压力低于液体的饱和蒸汽压时，进行液柱分离计算，即用下式判断液柱是否分离：

H_P,j+H_air-Z_j＜H_vapo_r(37)

其中，H_air为当地大气压力的液柱高度，单位为m；Z_j为节点j距基准面的高度，单位为m；H_vapor为液体饱和蒸汽压力，单位为m。

如图10所示，为液柱分离边界示意图。液柱分离边界的特点是该点压力相同，皆为油品的饱和蒸汽压，但是上下游流量不同，液柱分离边界节点的计算方法为：

其中，式中，分别为液柱分离点上游边界在j+1时刻的压头和流量； 分别为泄漏点下游边界在j+1时刻的压头和流量；分别为液柱分离点上、下游边界在j时刻的流量；H_vapor为液体饱和蒸汽压力，单位为m；分别为液柱分离点在j时刻的速度；ψ为加权因子，且0.5<ψ<1。若说明液柱依然分离，继续按液柱分离进行计算；若说明液柱已经合拢，此节点恢复为一般的内节点，按内节点计算。

7)根据步骤5)中得到的管内油流控制方程的离散格式以及步骤6)中确定的各边界条件对管道全线进行水热力耦合模拟计算，得到泄漏过程中所需相应计算时间内泄漏管道的泄漏量。

如图11所示，对管道全线进行水热力耦合模拟计算时，包括以下步骤：

①首先确定管道全线的基本参数，包括管道直径、管长、管道里程和高程。

②根据管道全线的基本参数，确定该管道物理模型中的相应物理单元，根据步骤6)中的边界条件处理方法，设置相应的边界条件和初始条件。

③根据当前温度，计算管道内油品的物性参数，包括油品密度ρ、粘度μ和列宾宗摩阻系数f。

④判断当前节点是否发生液柱分离，若发生液柱分离，则采用液柱分离的边界条件计算该节点当前时步的压力和流量值；若没有发生液柱分离，则根据当前时步各节点的油品物性，计算管道各节点当前时步的压力和流量值。

⑤根据各节点当前时步的压力和流量值，采用v特征方程计算各节点油品当前时步的温度值。

⑥判断所得到的温度值是否收敛，如果是，表明当前时步下各管段节点的温度均计算完毕，进入步骤⑦，如果否，则返回步骤③，将得到的温度值返回继续进行油品物性计算，直到收敛。

⑦确定当前时步油品批次界面位置，并根据得到的温度值以及泄漏点的边界条件计算当前时步的泄漏量。

⑧对以后每一时步重复步骤③～⑦，直至达到设定的计算时间。

下面结合具体实施例，对本发明方法进行进一步阐述。本实施例中，分别于2016年7月13日和8月13日在西南成品油管道“贵阳—安顺”段进行泄漏实验。

(1)实验过程

如图12所示，为本实施例中“贵阳-安顺”管道纵断面图。本实施例依托西南成品油管道“贵阳-安顺”管段改线工程，将贵阳与安顺之间的金银山改线点作为泄放点，实验管段全长93.889km，泄放点位于90.5km处。将油品泄放到油罐车中以模拟泄漏过程。泄漏实验通过管道打孔，并以软管将油品泄放至油罐车实现。具体实验结果如表1所示。

表1贵阳-安顺管道泄漏实验数据

	第一组实验	第二组实验	第三组实验
				泄漏点位置	90.5km	90.5km	84.2km
实验开始时间	13:26:00	12：10：29	15：41：40
				实验结束时间	13:45:15	12：26：07	15：59：47
阀门开度	15％	30％	30％
				油品泄漏量现场测定值	3.662m3	4.350m3	3.536m3

(2)计算结果对比

根据三次实验过程中上游的压力和流量数据以及下游的压力数据，对三次实验工况下的泄漏量进行测算，得到三次实验中的计算结果与实验结果对比如表2所示。

表2三次计算结果与实验记录数据对比分析

次数	现场泄漏量/m3	测算泄漏量/m3	阀门开度	相对误差(％)
					1	3.662	3.88	15％	5.9
2	4.35	4.94	30％	13.5
					3	3.53	3.38	30％	4.2

由表2可知，三次实验的误差都在15％以内，满足工程实际的计算精度要求。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (10)

1.一种液烃管道泄漏量测算方法，特征在于其包括以下步骤：

1)首先建立泄漏管道的通用物理模型；

2)根据建立的泄漏管道的通用物理模型，建立泄漏管道中各管段的连续数学模型，得到描述管内油流的控制方程；

3)采用各管段统一时步法对泄漏管道和时间域进行网格划分，得到离散的计算区域；

4)采用特征线法，将步骤2)中建立的泄漏管道的连续数学模型进行离散化，得到描述管内油流控制方程的特征方程；

5)根据步骤3)中的网格划分结果，将步骤4)中得到的管内油流控制方程的特征方程转换为离散格式；

6)根据步骤1)建立的通用物理模型中各管段的相邻处不同的边界类型，确定不同边界条件，并建立其离散格式；

7)根据步骤5)中得到的管内油流控制方程的离散格式以及步骤6)中确定的各边界条件对管道全线进行水热力耦合模拟计算，得到泄漏过程中所需相应计算时间内泄漏管道的泄漏量。

2.如权利要求1所述的一种液烃管道泄漏量测算方法，其特征在于：所述步骤2)中，管内油流的控制方程包括管内油流的连续性方程、动量方程和换热方程；其中，管内油流的连续性方程为：

$\frac{\&part;\left(\mathrm{\&rho;}A\right)}{\&part;t}+\frac{\&part;\left(\mathrm{\&rho;}Av\right)}{\&part;x}=0;$

管内油流的动量方程为：

$\frac{\&part;v}{\&part;t}+\frac{1}{\mathrm{\&rho;}}\frac{\&part;p}{\&part;x}+v\frac{\&part;v}{\&part;x}=-gsin\mathrm{\&theta;}-\frac{\mathrm{\&lambda;}}{2D}v\left|v\right|;$

管内油流的换热方程为：

$\frac{\&part;T}{\&part;t}+\frac{{\mathrm{Ta}}^2{\mathrm{\&alpha;}}_p}{c}\frac{\&part;v}{\&part;x}+v\frac{\&part;T}{\&part;x}=\frac{{\mathrm{\&lambda;v}}^3}{2Dc}-\frac{4q}{\mathrm{\&rho;}Dc};$

$a=max\left(a_{o,g}\right),a_{o,g}=\frac{1}{\sqrt{{\mathrm{\&rho;}}_o(\frac{1}{k_o}+\frac{D_g}{E_g{\mathrm{\&delta;}}_g})}};$

式中，t为时间，单位为s；ρ为油品在管截面上的平均密度，单位为kg/m³；x为距泄漏管段起点的距离，单位为m；v为管内油品的平均速度，单位为m/s；A为泄漏管段截面积，单位为m²；g为重力加速度，单位为m/s²；θ是泄漏管段与水平方向的夹角，单位为rad；D为管内径，单位为m；λ为达西摩阻系数；q为油品与单位面积管壁单位时间内的热流密度，单位为W/m²，T为管内油品平均温度，单位为℃；α_p为油品的膨胀系数，单位为℃^-1；c为油品的比热容，单位为J/(kg·℃)；a_o,g为压力波在不同介质中的传播速度，单位为m/s；ρ_o为不同批次的油品密度，单位为kg/m³；k_o为对应批次油品的弹性模量，单位为Pa；D_g为不同管段直径，单位为m；E_g为对应管段的杨氏弹性模量，单位为Pa；δ_g为对应管段的壁厚，单位为m。

3.如权利要求1所述的一种液烃管道泄漏量测算方法，其特征在于：所述步骤4)中，得到的管内油流控制方程的特征方程为：

C⁺特征方程：

$\left\{\begin{array}{c}\frac{dx}{dt}=v+a\\ \frac{dv}{dt}+\frac{1}{\mathrm{\&rho;}a}\frac{dp}{dt}=-gsin\mathrm{\&theta;}-\frac{\mathrm{\&lambda;}}{2D}v\left|v\right|\end{array}\right.;$

C^-特征方程：

$\left\{\begin{array}{c}\frac{dx}{dt}=v-a\\ \frac{dv}{dt}-\frac{1}{\mathrm{\&rho;}a}\frac{dp}{dt}=-gsin\mathrm{\&theta;}-\frac{\mathrm{\&lambda;}}{2D}v\left|v\right|\end{array}\right.;$

v特征方程：

$\left\{\begin{array}{c}\frac{dx}{dt}=v\\ c\frac{dT}{dt}-\frac{{\mathrm{T\&alpha;}}_p}{\mathrm{\&rho;}}\frac{dp}{dt}-\frac{{\mathrm{\&lambda;v}}^3}{2D}=-\frac{4K}{\mathrm{\&rho;}D}\left(T_0-T\right)\end{array}\right..$

4.如权利要求1所述的一种液烃管道泄漏量测算方法，其特征在于：所述步骤5)中，得到的所述控制方程的离散格式为：

C⁺特征方程：

$\left\{\begin{array}{c}\frac{\mathrm{\&Delta;}x}{\mathrm{\&Delta;}t}\&ap;+a\\ \frac{a}{gw}\left(Q_P^{j+1}-Q_A^j\right)+\left(H_P^{j+1}-H_A^j\right)+{\mathrm{fQ}}_P^{j+1}\left|Q_A^j\right|a\mathrm{\&Delta;}t=0\end{array}\right.;$

C^-特征方程：

$\left\{\begin{array}{c}\frac{\mathrm{\&Delta;}x}{\mathrm{\&Delta;}t}\&ap;-a\\ \frac{a}{gw}\left(Q_P^{j+1}-Q_B^j\right)+\left(H_P^{j+1}-H_B^j\right)+{\mathrm{fQ}}_P^{j+1}\left|Q_B^j\right|a\mathrm{\&Delta;}t=0\end{array}\right.;$

v特征方程：

$\left\{\begin{array}{c}\frac{\mathrm{\&Delta;}x}{\mathrm{\&Delta;}t}=v\\ c(T_P^{j+1}-T_A^j)-\frac{0.5(T_P^{j+1}+T_A^j)}{\mathrm{\&rho;}}{\mathrm{\&alpha;}}_P(P_P^{j+1}-P_A^j)-\frac{f}{{\mathrm{DA}}^3}{Q}_P^{j+1}{Q}_A^j(Q_P^{j+1}+Q_A^j)\&times;0.5\mathrm{\&Delta;}t+\\ \frac{4K}{\mathrm{\&rho;}D}\&lsqb;0.5\&times;(T_P^{j+1}+T_A^j)-T_0\&rsqb;\mathrm{\&Delta;}t=0\end{array}\right.;$

式中，Δx为空间步长，单位为m；Δt为时间步长，单位为s；f为列宾宗摩阻系数；w分别为管道横截面，单位为m²；分别为节点P处j+1时刻的流量，压头，温度和压力；分别为节点A处j时刻的流量，压头，温度和压力；分别为节点B处j时刻的流量和压头。

5.如权利要求1所述的一种液烃管道泄漏量测算方法，其特征在于：所述步骤6)中，得到的边界条件包括上下游边界、阀门边界、变径点边界、泄漏点边界、混油边界和液柱分离边界。

6.如权利要求5所述的一种液烃管道泄漏量测算方法，其特征在于：

所述上下游边界中，上、下游边界节点的压头和流量计算公式分别为：

$H_{i,N}^{j+1}=R_{i,N-1}^j-S_{i,N-1}^j{Q}_{i,N}^{j+1};$

$H_{i+1,0}^{j+1}=R_{i+1,1}^j+S_{i+1,1}^j{Q}_{i+1,0}^{j+1};$

式中，分别为管段i的第N个节点处j+1时刻的压头和流量，也即管段i上游边界的压头和流量；分别为管段i的第N-1个节点处j时刻的压头和流量方程的系数；分别为管段i+1的第0个节点处j+1时刻的压头和流量，也即管段i上游边界的压头和流量；分别为管段i+1的第1个节点处j时刻的压头和流量方程的系数；

所述阀门边界中，阀门边界节点的压头和流量的计算公式分别为：

$H_{i,N}^{j+1}=R_{i,N-1}^{+j}-S_{i,N-1}^{+j}{Q}_{i,N}^{j+1};$

$H_{i+1,0}^{j+1}=R_{i+1,1}^{-j}+S_{i+1,1}^{-j}{Q}_{i+1,0}^{j+1};$

$Q_{i,N}^{j+1}=Q_{i+1,0}^{j+1}=Q_P;$

$H_{i,N}^{j+1}-H_{i+1,0}^{j+1}={\mathrm{KQ}}_P^2;$

式中，分别为阀门上游边界在j+1时刻的压头和流量；分别为阀门上游边界j时刻的压头和流量方程的系数；分别为阀门下游边界j+1时刻的压头和流量；分别为阀门下游边界j时刻的压头和流量方程的系数；Q_p为阀门处的流量；

7.如权利要求5所述的一种液烃管道泄漏量测算方法，其特征在于：所述变径点边界中，变径点处压头和流量的计算公式为：

$H_{i,N}^{j+1}=R_{i,N-1}^{+j}-S_{i,N-1}^{+j}{Q}_{i,N}^{j+1};$

$H_{i+1,0}^{j+1}=R_{i+1,1}^{-j}+S_{i+1,1}^{-j}{Q}_{i+1,0}^{j+1};$

$H_{i,N}^{j+1}=H_{i+1,0}^{j+1}=H_P;$

$Q_{i+1,0}^{j+1}=Q_{i+1,0}^{j+1}=Q_P;$

式中，分别为变径点上游边界j+1时刻的压头和流量；分别为阀门上游边界j时刻的压头和流量方程的系数；分别为阀门下游边界j+1时刻的压头和流量；分别为阀门下游边界j时刻的压头和流量方程的系数；H_p、Q_p分别为变径点上下游处的压头和流量；

8.如权利要求5所述的一种液烃管道泄漏量测算方法，其特征在于：所述泄漏点边界中，泄漏点处压头和流量的计算公式为：

$H_{i,N}^{j+1}=R_{i,N-1}^{+j}-S_{i,N-1}^{+j}{Q}_{i,N}^{j+1};$

$H_{i+1,0}^{j+1}=R_{i+1,1}^{-j}+S_{i+1,1}^{-j}{Q}_{i+1,0}^{j+1};$

$H_{i,N}^{j+1}=H_{i+1,0}^{j+1}=H_P;$

$Q_{i+1,0}^{j+1}+Q_{i+1,0}^{j+1}+q_P=0;$

且泄漏点处压头H_p为：

$H_P=\frac{1}{A^2}\&lsqb;(AB+{\mathrm{gC}}_0^2)-\sqrt{(AB+{\mathrm{gC}}_0^2)-A^2(B^2+2{\mathrm{gC}}_0^2{H}_0)}\&rsqb;;$

式中， 分别为泄漏点上游边界j+1时刻的压头和流量；分别为泄漏点上游边界j时刻的压头和流量方程的系数；分别为泄漏点下游边界j+1时刻的压头和流量；分别为泄漏点下游边界j时刻的压头和流量方程的系数；H₀为泄漏孔外部的压头，单位为m；C₀为泄漏孔的流量系数，且C₀＝αC_dw≈0.6～0.65w；q_P为泄漏孔处瞬时泄漏流量，单位为m³/s；

9.如权利要求5所述的一种液烃管道泄漏量测算方法，其特征在于：所述混油边界中，混油边界的压力和流量计算公式为：

$P_{i,N}^{j+1}={\mathrm{\&rho;}}_{B}g(R_{i,N-1}^{+j}-S_{i,N-1}^{+j}{Q}_{i,N}^{j+1}-Z);$

$P_{i+1,0}^{j+1}={\mathrm{\&rho;}}_{A}g(R_{i+1,1}^{-j}+S_{i+1,1}^{-j}{Q}_{i+1,0}^{j+1}-Z);$

$P_{i,N}^{j+1}=P_{i+1,0}^{j+1}=P_P;$

$Q_{i+1,0}^{j+1}=Q_{i+1,0}^{j+1}=Q_P;$

$Q_P=\frac{{\mathrm{\&beta;R}}_{i,N-1}^{+j}-R_{i+1,1}^{-j}+(1-\mathrm{\&beta;})Z}{{\mathrm{\&beta;S}}_{i,N-1}^{+j}+S_{i+1,1}^{-j}};$

其中，ρ_A为前行油品的密度，kg/m³；ρ_B为后行油品的密度，kg/m³；β＝ρ_A/ρ_B；Z为混油边界处的高程，单位为m；分别为混油界面点上游边界j+1时刻的压头和流量；分别为混油界面点上游边界j时刻的压头和流量方程的系数；分别为混油界面点下游边界j+1时刻的压头和流量；分别为混油界面点下游边界j时刻的压头和流量方程的系数；P_p、Q_p分别为混油界面点的压力和流量；

所述液柱分离边界中，液柱分离边界节点的计算方法为：

$H_{i,N}^{j+1}=H_{i+1,0}^{j+1}=H_{vapor};$

$V_P^{j+1}=V_P^j+\&lsqb;\mathrm{\&psi;}(Q_{i,N}^j-Q_{i+1,0}^j)+(1-\mathrm{\&psi;})(Q_{i,N}^{j+1}-Q_{i+1,0}^{j+1})\&rsqb;\mathrm{\&Delta;}t;$

式中，分别为液柱分离点上游边界在j+1时刻的压头和流量；分别为泄漏点下游边界在j+1时刻的压头和流量；分别为液柱分离点上、下游边界在j时刻的流量；H_vapor为液体饱和蒸汽压力，单位为m；分别为液柱分离点在j时刻的速度；ψ为加权因子，且0.5<ψ<1。

10.如权利要求1所述的一种液烃管道泄漏量测算方法，其特征在于：所述步骤7)中，对管道全线进行水热力耦合模拟计算的方法包括以下步骤：

①首先确定管道全线的基本参数，包括管道直径、管长、管道里程和高程；

②根据管道全线的基本参数，确定该管道物理模型中的相应物理单元，根据步骤6)中的边界条件处理方法，设置相应的边界条件和初始条件；

③根据当前温度，计算管道内油品的物性参数，包括油品密度ρ、粘度μ和列宾宗摩阻系数f；

④判断当前节点是否发生液柱分离，若发生液柱分离，则采用液柱分离的边界条件计算该节点当前时步的压力和流量值；若没有发生液柱分离，则根据当前时步各节点的油品物性，计算管道各节点当前时步的压力和流量值；

⑤根据各节点当前时步的压力和流量值，采用v特征方程计算各节点油品当前时步的温度值；

⑥判断所得到的温度值是否收敛，如果是，表明当前时步下各管段节点的温度均计算完毕，进入步骤⑦，如果否，则返回步骤③，将得到的温度值返回继续进行油品物性计算，直到收敛；

⑦确定当前时步油品批次界面位置，并根据得到的温度值以及泄漏点的边界条件计算当前时步的泄漏量；

⑧对以后每一时步重复步骤③～⑦，直至达到设定的计算时间。