基于天然气放空系统的水力计算软件开发

娄 晨，林 棋

中国石油大学（北京）油气管道输送安全国家工程实验室 （北京 昌平 102249）

高压天然气放空系统由干线内待放空管存气体、放空管路和外界大气组成（图1），其中放空管路包括相应的管件及附件（包括管段、阀门、弯头、三通及稳固定支架等）。管段内气体流动属于典型的非定常流动，放空过程将经历超临界流、临界流及亚音速流3种状态，管段水力计算不可简单套用达西公式。相关的研究文献内容较少，目前关于天然气放空系统的放空量计算多采用一个没有依据的公式 （即管容、放空时间），计算结果与现场结果差异极大，导致工程设计与工程实际严重不符，此种情况的存在将带来严重的安全隐患 （比如，放空点火的热辐射伤害）[1-2]，故有必要开发出可适用于工程实际的放空管路水力计算软件。基于Fano方程，通过C语言编写开发出水力计算软件，并利用图解法与该软件进行瞬时放空量计算对比。同时结合东河伴生气管线下沉工程放空现场数据，进一步验证了计算软件的准确性。此软件的开发将为今后输气干线现场放空作业提供切合实际、安全、高效的放空方法，同时也为相关的科研研究提供便利的参考依据。

图1 高压天然气放空系统示意图

1 天然气放空管路水力系统理论

1.1 放空过程状态描述

高压天然气放空系统的实现是在短距离高压差的条件下，其放空过程非常剧烈，放空管段沿线的气体基本参数（温度T、密度ρ、速度V、压力p、马赫数Mach等）变化快、差异大，该过程为非定常流动，需要经历3个流动状态：超临界流（即：雍塞流）、临界流、亚音速流。

1）超临界流：放空初始时刻，放空干线内压力p0大，从而使得入口截面处压力p1很高，此时出口截面压力p2远大于外界大气压，放空管路出口处的气体流速达到当地声速 （即：临界流速，此时马赫数Mach2为1），对于给定的放空设施，此时将达到最大瞬时放空量（即：质量流量）。出口截面的气体依靠剩余压力差（p2-pa）进行膨胀，直至压力降至外界大气压。随着超临界流的时间推移，放空管路进出口压力均逐渐降低，但在出口压力p2降至大气压之前，此放空过程仍处于超临界流状态，在新的压力条件下放空过程达到新的最大瞬时放空量。

2）临界流：随着超临界流的进行，p0、p1及 p2均逐步减小，当出口截面压力p2下降到外界大气压力p0时，超临界流结束，此刻便为临界流状态，此时放空管路出口处的气体流速为临界流速（Mach2为1）。

3）亚音速流：当放空过程越过临界流继续泄压时，气体进入亚音速流状态。入口截面压力p1继续降低，出口截面压力p2保持不变（等于pa），出口气体流速逐渐减小（Mach2<1），质量流量也随之下降，在放空干线内压力p0下降至pa之前，该流动状态将持续进行，直至天然气放空系统内各压力值相等（pa=p0=p1=p2），放空过程结束。

1.2 天然气放空系统基本方程推导

设放空过程中在微元时间∂t内，气源滞止压力p0保持不变，在等截面水平放空管路中取∂x微元段作为控制体[3]。在此条件下的控制体流动基本方程（质量守恒方程、能量守恒方程及动量守恒方程）：

式中：m 为单位面积质量流率，kg/m2·s；ρ为气体密度，kg/m3；i为单位质量的气体热焓，J/kg；vi为放空系统不同截面的气体速度，m/s；d为放空管内径，m；f为水力摩阻系数；V为气体流速，m/s；式中下标1表示放空管路入口截面，下标2表示放空管路出口截面。

其他所需基础方程 （气体状态方程及相关热力方程）：

式中：Ti为放空系统不同截面的气体温度，K；γ为气体的等熵指数；pi为放空系统不同截面的气体压力，Pa；Mach为气体马赫数；M为气体摩尔质量，kg/mol；R为气体常数；a为放空系统所处的当地声速，m/s；Cp为气体的定压比热，J/kg·K；Z 为截面处气体压缩因子。

将方程组（2）代入方程组（1）中的质量守恒方程可得到关联方程组（3）：

同理将方程组（2）代入方程组（1）中的能量方程式可得关系式：

将关系式（4）代入关联方程组（3）可将其转化为马赫数与压力值的关系式：

同理将方程组（2）代入方程组（1）中的动量守恒方程可得关系式：

将方程组（3）及关系式（5）进行积分（考虑到气体等熵指数为常数），并将结果代入方程（6）可得：

由科尔布鲁克摩阻系数计算式（8）所示，管路中摩阻系数取决于雷诺数Re及管内壁粗糙度，在放空管路中，高速泄压的高流速使得管路雷诺数Re非常大，进而其对管路沿线的摩阻系数影响很小，在此为简化方程，忽略雷诺数Re的影响（即：管路摩阻系数仅由选定管路的内壁粗糙度决定）[4]：

式中：μ为放空管路内壁绝对粗糙度，m。

由此可对上述方程式（7）沿放空管路长度L进行积分，可得：

气体从输气干线通过圆形扩孔进入放空管路（p0T0⇒p1T1）入口截面的热力过程，可以视为一个典型的等熵热膨胀过程，故由热力学基本公式可得如下关系式[5]：

上述方程 （9）为有摩擦绝热一维流动方程式（即：Fano方程），将上述10个方程式（组）作为天然气放空系统水力求解计算关联方程组。

2 放空管路水力计算软件开发

2.1 程序计算流程简述

1）放空管路基础数据。方程式（9）中左边项 fL/d（简记为U）为阻力项，该项的求解可以分为3部分，如式（11）所示：

其中入口及出口的阻力项可以通过查阅资料获取，且两者占总阻力项的比例也较小。放空管段中L为其等效总长度（包括：弯头、三通、阀门的等效长度以及放空管段本身管长），再结合选取管段的管内壁绝对粗糙度，即可求得Fano方程阻力项。

2）超临界流初始时刻质量流量（瞬时放空量）。将方程式（7）中的气体等熵指数γ视为定值（Cp/Cv=1.3），放空初始时刻放空管路出口处的气体流速为临界流速（Mach2为 1），结合方程式（11）可以将方程式（7）转化为仅有Mach1的迭代计算式：

通过程序的迭代计算可以求得入口截面气体马赫数 Mach1，将其代入方程式（5）以及方程组（10）中的压力关系式，即可求得放空管路进出口压力p1、p2。将p1及Mach1代入方程组（10）中的温度关系式，可求得入口截面气体温度T1，再将T1、Mach1及Mach2代入方程式（4）中，可求得出口截面气体温度T2。经过上述求解后，可由式（13）计算求得超临界流的初始质量流量：

式中：W为瞬时放空量，kg/s；Z1为截面处天然气压缩因子。

3）临界流时刻质量流量（瞬时放空量）。由上述放空过程描述可知，当放空气体处于临界流时，出口截面压力p2始终等于大气压力pa，出口气体马赫数Mach2为1，由迭代计算式（12）可知，入口截面气体马赫数Mach1求解值不变，代入方程式（5）可得入口截面压力p1，再代入方程组（10）压力关系式，可求得临界流状态时干线管存压力值p0，由方程组（10）温度关系式可知入口截面气体温度求解值不变。将Mach1、Mach2及T1代入方程式（4）可求得出口截面气体温度T2，再同理利用方程式（13）可求得临界流时刻质量流量。

4）亚音速流质量流量（瞬时放空量）。当放空过程处于亚音速流时，出口截面压力p2保持不变（等于pa）。针对临界流状态计算所得的干线管存压力值，采用压力等分递减求解计算方法。取一个合适的压力递减值Δp0，由此可得下一时刻的p0，代入方程组（10）压力关系式，求得p1，由于此时迭代计算式（12）存在Mach1及Mach2两个未知量，故先假设Mach2=1，通过迭代计算求得Mach1，之后将Mach1作为已知值，将方程式（9）转化为求解Mach2的迭代计算式进行迭代求解，计算出Mach2。重复上述过程，直至上述两个迭代式所求得的Mach1及Mach2均达到稳定收敛状态，将最终时刻迭代值赋予Mach1及Mach2。此时同理由方程组 （10）温度关系式求得入口截面气体温度、由方程式（4）可求得出口截面气体温度T2、有方程式（13）可求得亚音速流任一压力递减值所对应的质量流量。

5）放空过程累计放空时间。对于天然气放空系统，最为关注的因素就是其完成整个放空过程的累计放空时间。上述介绍了各个流态下的质量流量计算思路，但仅由上述方法及关联方程组是无法求得累计放空时间，在此引入一个经典基本数值积分法——“梯形法”。将初始干线管存压力以某一微小压力递减值进行等分：

由2）～4）计算推导过程，可对上述每一个压力等分递减值计算出其所对应的瞬时放空量Wi，再结合方程组（2）及放空干线标准管存体积计算式，通过式（15）求得各个计算时刻的管存储气量Qi，由于微小压力递减值非常小，在每个计算区间可以看成匀速放空过程，由此可以求得各个压力区间段的等效放空时间ti，将各区间的放空时间求和即可求得放空过程的累计放空时间T。

式中：T为累计放空时间，h；Q为输气干线管段储气量，m3。

6）放空管路管径选取。上述推导了放空系统的水力计算详细计算步骤，由规定的总累计放空时间可以反推计算求得所需要选取的放空管路管径，由此可更好地运用在现场工程实际。反算流程如图2所示。

图2 部分C语言程序计算流程图

2.2 算例计算分析

为检验程序的准确性，利用所开发的软件 （图3），建立放空系统算例模型（表1），并同时采用图解法进行求解，API 521提供了以可压缩流体有摩擦绝热一维流动的Fano方程为基础的图解近似计算法（Lapple图）[6]。将两者结果进行对比验证。

图3 软件计算界面图

软件计算可求得放空系统各时刻气源压力值所对应的瞬时放空量、累计放空时间及单位质量流率、进出口截面的压力及温度值、各时刻管段储气量部分数据提取见表2。从中获知：

1）累计放空时间随管存压力值的降低逐渐增大，放空时间的增长速率也逐渐增大，且在亚音速流区域呈现出急速增长。

2）瞬时放空量及单位质量流率在泄压放空过程中呈现出线性递减的变化趋势，与放空过程所经历的3种流态形式没有关系。

3）放空管路入口截面压力p1及输气干线管存储气量随干线压力的减少呈现出与累计放空时间相反的变化规律，在雍塞流区域急速下降，而在亚音速流区域下降幅度平缓。

4）放空管路出口截面压力p2在雍塞流区域急速下降至大气压力pa，然后在亚音速流区域维持该压力状态不变。

5）放空管路入口截面温度T1略小于输气干线气体温度T0（由于开孔处微弱的节流），在雍塞流区域维持温度不变，在亚音速流区域缓慢升温至干线气体温度值，从而达到平衡（图4）。

表1 放空系统程序计算算例基本参数

6）放空管路出口截面温度T2受初始放空时刻的强烈节流效应，温度急速下降至某一最低温度，之后随放空过程的进行，逐渐升温，最终达到T2=T1=T0；此升温过程在雍塞流区域上升幅度较慢，而在亚音速流区域呈现出急速回温的变化趋势（图4）。

表2 放空系统部分计算结果汇总表

图4 放空系统内各温度变化曲线

7）在雍塞流状态下压降速度快，压降变化值占到总放空过程的90%左右，且持续时间较短（小于50%），瞬时放空量大，会使得放空管路出口处急剧降温（初始瞬间可达-40℃左右）；亚音速流状态下压降速度缓慢，压降变化值仅占到总放空过程的10%左右，但持续时间却大于雍塞流状态，放空完成时刻系统各压力值与大气压pa一致。

8）同时利用图解法对算例进行了求解，发现程序计算瞬时放空量与图解法计算结果误差基本控制在2%以内（图5）。

2.3 东河伴生气管线下沉工程放空现场

2013年9月底在南疆铁路二线建设中，轮台县至库车县A标志段B涵洞处，因通过涵洞下的公路路面基线与铁路桥距离太小，当地车辆通行不便，需要降低路面，使此处东河伴生气管线埋深不符合设计要求，因此需要进行伴生气管线沉降施工。需将东河伴生气管线出站至3号阀室进行放空作业。本工程的施工程序如下：原管线开挖—管线放空—氮气置换—新管线预制—焊缝无损检测—割管—封堵—新管线安装—焊缝无损检测—补口—电火花检漏—补伤—管线吹扫—氮气置换—管线复产—管线回填。

图5 图解法与计算软件结果对比

东河气站放空管线主阀为一电磁阀，通过系统操作全开时开度也不能达到100%，为提高放空速度，经协调通过打开旁通阀进行放空。放空管线长度34.4km，管径168×6/7mm。放空前，管线压力为6.15MPa，管段内天然气量标准状况下为40 408m3。放空管线规格DN100，放空阀门开度100%。在管线放空过程中，现场作业中工作人员对放空时间做了详细记录，相关数据见表3。

利用开发的水力计算软件，根据现场进行天然气放空系统建模，将程序计算结果与现场放空记录数据进行对比（图6），进一步验证本计算软件的准确性。

由对比结果可知：两者数据对比基本一致，放空过程压降变化曲线重合度较高，软件计算所得压降曲线较为圆滑规整，其超临界流持续时间稍短于实际放空现场，软件计算瞬时放空量稍大；现场记录数据显示，当旁通处于关闭时，放空速度较慢，打开旁通后，放空速度加快，压降变化曲线与软件计算基本一致，软件计算累计放空时间略小于实际放空现场。上述存在的差异主要为：①软件计算模型暂未考虑放空系统的等效长度（阀门、弯头及三通等）；②现场放空时，旁通打开的时间是在放空进行了一个多小时之后。由于这两点的差异性，导致了软件计算天然气放空速度略快。但从总体上看，软件计算精度较高，可以应用到实际现场。

表3 东河伴生气管线放空现场数据记录表

图6 放空现场与计算软件结果对比

3 结 论

1）基于C语言编写的天然气放空系统水力计算软件，操作方便、计算速度快，利于工程实际应用。

2）由软件计算算例分析可得：超临界流状态下压降速度快，压降变化值大，持续时间较短，瞬时放空量大，使得放空管路出口处急剧降温；亚音速流状态下压降速度缓慢,压降变化值小,但持续时间较长。

3）计算软件与图解法计算结果进行对比，精度较高（误差在2%以内），并弥补了图解法仅可求得瞬时放空量而无法得到其他参数的缺陷。

4）通过东河伴生气管线现场天然气放空所记录的数据，利用软件建模计算，对比结果显示软件计算结果与实际现场放空吻合度高，可将软件应用于工程实际，为现场放空作业提供准确、安全、高效的放空方案。

[1]周雪漪.计算水力学[M].北京：清华大学出版社,1995.

[2]API R P521 Guide for pressure relieving and de-pressuring system[S].

[3]叶学礼.天然气放空管路水力计算[J].天然气工业，1999,19（3）:90-94.

[4]李玉星,姚光镇.输气管道设计与管理[M].东营：中国石油大学出版社,2009.

[5]叶学礼.图解法求天然气瞬时放空量[J].天然气与石油,1999,17(2):1-3.

[6]K.K.BOTROS,W.M.JUNGOWSKI,M.H.WEISS.Models and MethodsofSimulating GasPipeline Blowdown[J].The CANADIAN JOURNAL OF CHEMICAL ENERGINEERING,1989，67（8）:529-539.