燃气管网置换过程动态仿真计算与研究1)

王安然 张 凯 ,2) 赵 瑶 苗 芊 杨荣超 曾 波 张 勍 史占东

*(中国计量大学，杭州 310018)

†(中国烟草标准化研究中心，郑州 450001)

目前燃气管网的投产置换大都采用无清管器氮气隔离的方式[1]，其主要分为两步：氮气置换和燃气置换[2-3]。在置换开始之前，通常会采用经验公式对置换过程所需的燃气量进行估算[4]，但在实际投产置换过程中，管网运行的工况不断发生变化[5]，而经验公式无法适应多变的工况，这会导致估算的燃气量与实际需要相差较大，此外，为确保安全和燃气浓度的达标，往往会持续输入大量燃气，过量的燃气会通过放空或燃烧的方式进行处理，导致大量的资源浪费，造成不必要的经济损失[6-7]。

对燃气管网置换过程中流体状态参数进行仿真计算是个较为复杂的研究，大多数研究讨论的是流体动力学方程的选取和求解，或是对单直管道进行置换过程的仿真计算。Johnson等[8]提出了能够用于计算气体置换过程中混合气体区域长度的计算公式，它可以应用于直接置换或惰性气体置换。郭谊民等[9]研究了气体在管道内流动的二维湍流扩散理论，研究了流速、置换时间与混气段长度的关系。张鹏云[10]采用实验的方法模拟输气管道氮气置换过程，得出了体积浓度与时间和长度的变化关系。段威等[11]推导出混气段长度公式，采用Fluent软件分析了介质流速、管径、管线长度对混气段长度的影响，但模型简化较多，模型结构只包含直管和弯管。Ma等[12]建立了一个理想化模型和一个近似于实际工况的模型，研究了如流体流入方式、末端尺寸及末端压力等不同条件对置换的影响并对结果进行对比分析。南军虎等[13]对流体流动中阻力损失产生的原因进行了清晰解释，对如何选取流体动力学方程进行理论计算作出了明确阐述。叶恒等[14]通过仿真计算得到中俄东线输气管道置换最优注氮管容比。Jeroen等[15]采用不同阶数的欧拉方程对管网进行仿真，讨论了不同模型的误差，并对模型的适应性进行了阐述。

算法研究方面，宋杨等[16]给出了能够计算没有环状支路管网的动态仿真算法，其算法具有计算敏捷、容易收敛的特点，但研究只停留在理论层面的推导证明，并未体现实际工程中的应用效果。Desideri等[17]提供了一种能够通过输入和输出条件（由压力流量相关性定义）来解决管网问题的算法，该方法特别适合求解边界条件由压力流量二次曲线定义的管网。近年来，管网算法的选取热衷于采用遗传算法对管网进行求解。张喜明等[18]采用遗传算法对燃气管网中的参数进行仿真计算，实现了对管径的优化设计并用来节省建设费用。高建丰等[19]使用改进的遗传算法大大加快了计算的收敛速度，但局限性在于考虑的因素较少，与实际工况存在一定差距。

现阶段，由于传统三维仿真建模难度大、模型灵活性差、计算时间长等特点，实际工程应用中大都采用一维管道网络模型对燃气管网进行仿真计算。但一维模型存在无法完全还原管道中的三维阻力结构的局限性，其精度往往偏低。本文在将燃气管网抽象成一维模型的基础上进行动态求解，增加流阻元件并根据实用流体阻力手册精确计算管网中的阻力损失。将计算结果与Fluent仿真结果进行对比，验证程序仿真计算结果的可靠性，并分析不同管长、置换速度条件下水平直管道终点混气浓度的变化情况。最后将计算模型应用于实际投产置换的案例中，计算置换过程中各放散口燃气气头的到达时间以及置换合格的时间。

1 数学模型

1.1 燃气管网一维仿真模型

目前关于燃气管网置换问题的仿真大都通过三维仿真软件进行，但对于复杂燃气管网的仿真往往存在建模复杂、网格数庞大、计算周期长等问题，且利用现有的仿真软件解决问题时模型的灵活性较差，若管网结构后续因运营维护发生变化，则原有三维模型将无法支持变化后实际结构的仿真计算。为解决上述问题，编写燃气管网一维仿真程序，在保证计算精度的情况下，极大程度地减少参与计算的方程数量，同时由于其可编程性，相比于现有仿真软件的灵活性较大，能够适应任何结构的燃气管网计算问题。在燃气管网的一维模型研究中，通常将复杂的三维管网抽象成一维的管道与节点的结构，即燃气管网是由若干管道和节点构成的，其基本模型结构如图1所示。

图1 燃气管网一维结构模型Fig.1 One-dimensional model of gas pipeline network

在燃气管网实际运行和置换过程中，管网整体温度变化较小，由温度变化引起的参数变化可忽略不计[20]，输送过程中没有发生相变，所以可认为是单相可压缩流体。为构建本模型而做出的假设可以表述如下[21]。

（1）燃气管网流体网络内部不存在动力源。

（2）节点内流体状态均匀（内部压力处处相等）。

（3）流动阻力仅考虑局部阻力与沿程阻力，并将其附加在管道阻力上。

（4）同一支路管道截面积不变，工质参数使用相连节点参数加权平均来表示。

特别地，节点又可以进一步划分为内部节点与外部节点。其中内部节点将参与程序的迭代计算，外部节点一般为供气端和末端用户节点，用来设置边界条件。

1.2 燃气物性的计算

根据实际燃气管网中燃气组成的检测结果，燃气组分中不同气体的摩尔分数为：CO(23.69%)，H2(22.75%)， CH4(2.27%)， CO2(10.31%)， N2(36.16%)，H2O(4.82%)。为了方便计算，将燃气视为单一气体，其分子量可以通过各组分混合后的平均分子量进行表示，计算公式为

式中，M为气体平均分子量，Mn为各组分的分子量，Cn为各组分的摩尔分数。

在程序迭代过程中，考虑气体的可压缩性，可以根据理想气体状态方程对节点上气体的密度进行迭代更新。混合气体的气体常数Rg同样可由各组分的气体常数Rn求得，计算公式为

通过计算，可以确定此种燃气的平均分子量为22.98 g/mol，平均气体常数为 360.52 J/(kg·K)。

1.3 燃气管网流体网络计算

燃气管网一维流场计算公式见式（3）～式（5），其中节点用i表示 (i= 1, 2, ··,n)，管道由其两端节点i和j进行表示（j= 1, 2, ··,n；且j≠i）。

式中，ρi为节点i的密度，Vi为节点i的容积，Rij为管道ij的阻力系数，Gij为管道ij的质量流量，λij为管道ij的摩擦阻力系数，lij为管道ij的长度，dij为管道ij的直径，ξij为管道ij包含的局部阻力系数，ρij为管道ij内部流体的平均密度，Aij为管道ij的截面积，Uij为管道ij的平均流速，Pi为i点的压力，Pj为j点的压力，Hij为节点i，j之间的宏观动能、势能及动力源产生的压力，t为时间，Dij为节点i，j之间的连通关系，具体数值为

1.4 燃气管网阻力系数的计算

对于燃气管网的流场计算，其计算的准确与否主要取决于管网各项阻力系数是否准确。对于整个管网的阻力系数，主要有沿程阻力系数和局部阻力系数两种，而局部阻力系数中，管道分支处的排出三通阻力系数的确定则是长久以来管网仿真的一大难题。为精确计算燃气管网中的三通阻力系数，程序根据实用流体阻力手册中提供的三通阻力系数表，通过二维线性插值算法将离散的三通阻力系数表变为适应多种工况的阻力系数数据集，相对于传统的查表方式或根据经验确定阻力系数的方式，计算精度得到了大幅提升。三通阻力系数表见表1，二维线性插值算法示意图见图2。

表1 排出三通管的三通阻力系数Table 1 Resistance coefficient of discharge tee

由图2可以列出二维线性插值算法的公式为

图2 二维线性插值示意图Fig.2 Schematic diagram of two-dimensional linear interpolation

式中，T和p分别为点的横坐标值和纵坐标值。

在实际工程应用中，由于管道老化、堵塞等原因，其阻力会发生变化。传统仿真软件很难对不同管道进行阻力系数的修正。而自编程序具有很强的灵活性，可以根据管网运行的实际数据，针对性地对每条管道的阻力系数进行修正，并将其嵌入程序计算中，进一步提高仿真计算的精度。

1.5 管道终点混气浓度计算方法

因燃气与氮气密度接近，管道径向浓度梯度很小，主要的气体传质过程都发生在沿管道轴向，因此可用一维对流扩散方程描述管道内断面的摩尔分布，其公式为

其中，K1为燃气的摩尔分数；v为置换速度，m/s；DT为气体扩散系数，m2/s；x为任意截面至起始接触面的距离，m；t为置换经过的时间，s。

为方便求解式(7)，引入变量S，其中L为管道全长。将其代入并求解方程可以得到置换过程燃气浓度与时间距离(S)的函数为

根据边界条件可得：t=L/v时，S= 0，K1= 0.5；t= 0时，管道内未通入燃气，因此K1= 0，S为无穷大。引入概率积分函数φ(S)=代入边界条件解得C2=0.5。将C1和C2代入方程可得管线终点燃气浓度K1与管内残留氮气浓度K2分别为

对流扩散系数可以通过Levenspiel[22]给出的关系图来确定，其描述了对流扩散系数DT与雷诺数Re及施密特数Sc的关系如图3所示。

图3 管道中的对流扩散系数Fig.3 Convection diffusion coefficient in pipe

结合燃气管网的流场计算，终点混气浓度程序计算流程如图4所示。

图4 混气浓度计算流程图Fig.4 Gas concentration calculation flow chart

2 计算结果分析

2.1 水平直管道计算结果

2 .1.1 水平直管道燃气置换程序计算及仿真验证

结合程序计算和Fluent软件对水平直管道燃气置换进行模拟仿真。管道模型长为1500 m，管道直径为0.6 m，绝对粗糙度为0.08 mm，燃气从入口处以1 m/s的速度充入，出口设置为压力出口，初始情况下管道内充满氮气。标准k–ε模型由于其适用范围广、计算准确性高等优点而得到广泛应用，置换过程中管道长度较长，管内气体流动状态通常为湍流，因此适合采用标准k–ε湍流模型进行计算，近壁区选用壁面函数法。为检验Fluent仿真的可靠性，首先对算例进行网格无关性验证，网格数量与压损计算结果见表2。

分析表2可知，当网格数量达到1 860 000后，若继续对模型网格进行加密，对流场计算精度的提升有限，且会消耗大量计算资源，因此采用Mesh 4进行仿真即可满足精度需求。程序计算与Fluent仿真的结果对比如图5所示。

表2 流场参数与网格数量的关系Table 2 Relationship between flow field parameters and number of elements

图5表明，通过程序计算得出的水平直管终点燃气浓度变化与Fluent仿真得到的浓度变化具有较好的一致性，验证了程序计算的可靠性。在置换过程模拟中(取燃气浓度1%～99%)，平均相对误差为4.9%，满足精度要求。

图5 水平直管道终点处燃气浓度仿真对比Fig.5 Comparison of simulation results of terminal concentration in horizontal pipeline

2 .1.2 不同质量流量下混气段发展过程和终点混气浓度

为了研究置换速度对混气规律的影响并验证模型的可靠性，对水平单直管进行置换过程的仿真计算。管道长度选择为1 500 m，管道直径为0.6 m，分别对流速为 0.5 m/s，0.75 m/s，1.0 m/s，1.25 m/s，1.5 m/s的情况进行程序仿真计算，计算结果如图6所示。

对置换过程中燃气流量进行积分，得到不同置换速度V下置换完成所需要的燃气量M如表3所示。

由图6可知，通入燃气的流速越快，置换速度越快，终点燃气浓度的变化速率越快。由对流扩散理论可知，置换速度越快，湍流扩散系数会越大，分子扩散得越快，混气段长度增加得越快。由表3可知，置换完成时所需的燃气总量随着置换速度的加快而略有减少，因此在实际置换过程中，为了置换过程安全快速进行，也为了减少燃气的损耗，可适当加快置换速度。

图6 不同置换速度下终点燃气浓度随时间的变化Fig.6 Variation of terminal concentration with time under different displacement rates

表3 不同置换速度下置换完成所需燃气量Table 3 Required gas mass at different displacement rates

2 .1.3 不同管道长度下混气段发展过程和终点混气浓度

在实际投产中，燃气管网可能长达几百甚至上千米，因此长度对混气规律的影响至关重要。在1 m/s的置换速度下，分别对管长为500 m，1 000 m，1 500 m，2 000 m，2 500 m，管道直径为0.6 m的管道进行程序仿真计算，探究不同管长下终点混气浓度的变化规律，计算结果如图7所示。

由图7可知，管道长度越长，管道终点燃气浓度的变化速率越慢，特别是每条曲线拐点处的斜率随管道长度的增加而减小。因为在置换速度相同的情况下，管道越长，管道内的混气段越长，起始接触面与气头的距离拉长，混气段的置换时间也随之增加，导致终点处的浓度上升速率减小。因此在实际置换中，对于管网中长度远、管道粗的管段，可以在其管段中间部分增设阀门平台分段置换，保证管网快速、安全地进行置换。

图7 不同管道长度下终点燃气浓度随时间的变化Fig.7 Variation of terminal concentration with time under different pipe lengths

综上，采用程序计算模型对燃气管网投产置换过程中终点燃气浓度进行仿真计算是可行的，由于模型的构建基于流体力学守恒方程，在算法上具有较强的通用性，因此能够适用于大部分常见燃气管网结构。此外，程序计算相较Fluent仿真计算机内存占用率低，且简化了建模过程和网格划分过程，因此对使用人员要求低，在实际工程应用中具有更好的普及性。

2.2 燃气管网三线联合置换

2 .2.1 管网基本结构

某地区为推动当地陶瓷企业产业升级、稳步发展，有效缓解当地陶瓷产业环保压力，特建设大型燃气管网用于集中供气。为了管网系统的安全运行，在实际投产之前必须对管网整体进行气体置换。待置换的燃气管线包含线路1，线路2和线路5三条，其沿线长度分别约为8 000 m，2 400 m和8 900 m，管道容积分别约为24 300 m3，5 800 m3和23 250 m3。将其抽象之后的拓扑结构模型如图8所示。因为三条线路包含并线的管段，为了提高置换效率，采用三条线路联合置换的方式进行置换。置换开始前管道内充满氮气，置换起点为节点1，三条线路的置换终点均设有两个放散口，各节点的标号见图8。

图8中，节点1为置换起点，节点7，8，14，15，21，22为三条线路置换终点处的放散口，节点4～8是属于线路2的节点结构，节点16～22是属于线路5的节点结构，节点10～15是属于线路1的节点结构，剩余节点均处于并线的管线结构中。2F05和2F04对应线路2末端的放散口，5F14和5F11对应线路5末端的放散口，1F03和1F05对应线路1末端的放散口，1F02和5F05分别为线路1和线路5沿线途中的放散口。各管段的内径D和管长L的参数如表4所示。

表4 燃气管网结构参数Table 4 Structural parameters of gas pipeline network

图8 燃气管网结构图Fig.8 Structure diagram of gas pipeline network

2 .2.2 三线合并置换方案

首先对最近的线路2进行置换，即对3～4，4～5，5～6所在管线进行置换，6～7，6～8为阀门平台上的两个放散口，相对高度为5 m。当线路2置换一段时间之后，同时对9～10，10～11，11～13所在的线路1进行置换，11～12，13～14，13～15分别是管线沿途和管线末端的三个放散口，对应编号1F01，1F03和1F05，相对高度为5 m。当线路1进行置换一段时间后，同时对9～16，16～17，17～19，19～20所在的线路5进行置换，17～18，20～21，20～22分别是线路5沿途和末端的三个放散口，分别对应编号5F05，5F04和5F11，相对高度为5 m。为提高置换效率，加快混气段的置换速度，当气头到达1F02时，等待一段时间之后再关闭1F02。同理，当气头到达5F05后，待置换继续进行一段时间后再进行关闭。对于三条线路末端的6个放散口，需要分别待其浓度达标后（燃气含量95%以上）立即关闭。实际置换过程中的放散口如图9所示。

图9 燃气置换过程中的放散口Fig.9 Vent during gas replacement

2 .2.3 置换时间计算结果及实际置换时间对比

由于此实际管网用于对多家用户企业的集中供气，因此在厂内设置有大型放空火炬装置，用于将生产的多余燃气进行放空和燃烧处理。得益于此结构的存在，使得厂外用于供气的燃气管网的压力可以保持相对稳定，又因为在置换时由于整个置换过程是阶段性的，所以置换速度也会随之进行变化，为了更加真实准确地模拟实际置换情况，在供气压力维持在50 kPa左右的前提下，边界采用速度入口条件进行仿真计算，燃气管网内部各个节点的压力将会通过程序瞬态计算得出。其中，供气端的流速在不同时间节点的变化情况如表5所示。

表5 不同时间节点对应管道流速Table 5 Velocity of pipeline at different time

基于以上模型，仿真计算得到燃气气头到达各放散口的时间以及各放散口燃气含量达到95%的时间，并以该时间作为置换合格的时间，将仿真计算得到的气头到达时间和置换完成时间与现场实测数据进行对比，如表6所示。

表6 实际置换时间与程序计算置换时间Table 6 Actual displacement time vs.simulated displacement time

为了更加清楚地描述仿真计算与现场实测数据的偏差，对相对误差定义为

由相对误差计算公式可以得到各个放散口处气头到达时间与置换合格时间的相对误差，如表7所示。

表7 程序计算结果相对误差Table 7 Relative error of program calculation results

由表7可知，仿真计算结果与现场数据之间的误差范围在7.9%～13.7%之间，燃气气头到达时间和置换合格时间的平均相对误差分别为10.8%和11.6%，此误差可以满足工程上的需求，因此说明采用本系统进行燃气管网置换过程的动态仿真计算并进行混气分析与方案制定的方法是可行的。

对置换过程中燃气浓度在管网中的轴向分布进行分析时，选取线路1直管段较长的2～3，3～9，9～10，10～11管段进行分析。由表可知燃气气头到达时间为置换开始后1 860 s，途中沿线放散口1F02关闭的时间为3 360 s，置换期间10～11管段的平均流速为2.227 m/s，置换进行至620 s，1 240 s，1 860 s和3 360 s时线路1主管内的燃气摩尔分数分布规律如图10所示。

由图10可知，当置换进行至620 s时，混气段长度约为2 300 m；当置换进行至1 240 s时，混气段长度约为3 300 m；当置换进行至1 860 s时，燃气气头刚刚到达1F02放散口，此时管内混气段长度约为4 400 m；保持1F02打开继续放散至3 360 s时，混气长度约为3 700 m。在图10(d)中，管线内的燃气摩尔分数在混气段长度为7 600 m后发生陡降，其原因是混气段经过10～11管段时，由于1F02放散口的口径较大，大部分混合气体经放散口被排出主管，因此1F02放散口后的燃气浓度较低，使得混气段长度在此处得到抑制。

图10 不同时刻下线路1管道内的燃气浓度分布Fig.10 Gas concentration distribution of line1 at different times

由表6可知，8:18分即线路1置换开始3 360 s后，1F02被关闭，此时10～11管段的流速下降为0.607 m/s，线路1中的燃气浓度增加速率降低。节点11处的浓度随时间变化的曲线如图11所示。

图11 节点11处燃气浓度随时间变化曲线Fig.11 Variation curve of gas concentration with time at node 11

综上所述，在置换过程中，主管内的混气段长度会随着置换的进行而增加，但是由于沿途放散口的存在使得主管内的混合气从放散口排出，从而有效抑制混气段的长度。放散口关闭后，会导致置换速度下降，从而拖慢整个置换过程。根据以上数据与分析可归纳出此类置换工程的一些指导方案。

（1）在置换开始时，为加快置换进度，应尽可能多地打开放散口进行置换。

（2）若线路中前段存在大口径放散口，会有效抑制混气段的长度，但会拖慢后续管段的置换进度，同时考虑到现场人员安全以及污染物排放问题，应在气头到达时及时将其关闭。若线路中后段存在大口径放散口，应在燃气分数达到一定程度时再关闭，加快置换进程。

（3）在长度较长的管段部分如管道10～11，可以适当在中间设置放散口，以有效缩短混气长度。

（4）在工作人员安排上，由于线路2置换完成的时间与其他两线路置换完成时间的间隔较大，可在线路2置换完成后将其工作人员分配去其他未置换完成的线路，最大程度节约人力资源。

3 结论

本模型将一维管网仿真与流阻元件建模相结合，提出了一种新的燃气管网置换过程动态仿真计算方法。在保留一维模型计算速度快、建模简单等优势的基础上，采用二维线性插值算法得到准确的局部阻力系数，提高了模型的计算精度。由于模型计算利用的是管网参数的守恒关系，并将燃气管网中的结构进行了抽象和等效，因此能够适用于大多数燃气管网的动态仿真研究工作。

本文首先用水平直管示例验证模型的正确性和准确性，并得到置换过程中管网参数的变化规律：管道越长，终点混气浓度变化梯度越平缓，由于管道体积的增加，置换完成时所需通入的燃气量也会随之增加；若管道长度不变，流速越快，终点混气浓度变化梯度越陡，置换完成时所需燃气总量也会随之减少。在实际燃气管网置换工程的应用中，随着置换时间的增加，管网内的混气段长度会随之增加，使用沿线途中的放散口进行放散可以有效抑制混气段长度的增加，同时提高置换效率。将计算结果与实测数据进行对比，得到燃气气头到达时间的平均相对误差为10.8%，置换合格时间的平均相对误差为11.6%，满足实际工程需求，从而能为燃气管网置换方案的制定提供重要指导依据。