温凯，焦健丰，袁运栋，董楠，殷雄，樊迪，宫敬

1.中国石油大学（北京）城市油气输配技术北京市重点实验室；2.中油国际管道有限公司；3.中国石油工程建设有限公司

0 引言

随着“双碳”目标带动的能源结构转型持续推进，天然气作为最清洁的化石资源，在中国实现“碳达峰、碳中和”任务目标过程中起到桥梁和支撑作用［1-3］，其在中国能源消费占比从 2010年的 4%，提升至2020年的8.2%［4-5］。国家石油天然气管网集团有限公司（简称国家管网集团）的成立标志着整个天然气市场向着上游多主体供应、中游单管网集输、下游多市场竞争的“X+1+X”格局加速推进［6-7］。随着天然气在能源消费中的占比日渐提高、交易环节和主体的增多［8-10］，输差作为天然气管输过程中一个极其重要的指标越发引起重视。输差与天然气各环节主体的经济效益直接挂钩，受多种实际因素影响。为提高天然气贸易的精细化管理，保障天然气管网高效运行，有必要研究天然气管道输送的输差问题。天然气管道输差的产生原因分为计量偏差和计算不精确两大类［11-13］。近年来，研究人员在传统分析方法的基础上，借助人工智能、大数据分析等新算法深入分析输差的影响因素，最终，都将问题的关键点指向管存计算［14-19］。

目前，管存的计算方法主要有理论计算、软件仿真及在线数据分析等。研究人员基于实际气体的状态方程，采用理论公式计算管存的过程中发现公式法耗时费力，同时存在压缩因子计算不准确的问题［20-27］。林棋等［28］采用微粒群进阶算法改进了上述问题，同时制订了在不同输量下的管存控制原则。林杨等［29］在西气东输天然气管道的管存计算中设计了SCADA（数据采集与监视控制）数据质量补救方法，有效提高了计算准确度，并使用苏霍夫公式结合相关数据，得出实际的沿线温度，一定程度上提高了数据的可靠性。上述管存计算均基于稳态算法实现，但在管输过程中天然气往往处于非稳态。李猷嘉［30］比较了稳态与非稳态两种计算方法，结果表明稳态计算方法的结果偏于安全。

自20世纪90年代以来，国外各水力仿真模拟软件的引入为管存计算带来了帮助。邢晓凯等［31］将TGNET（气体管道瞬态和稳态模拟计算软件）模拟的管存值同稳态计算公式进行对比，结果表明稳态计算公式有所偏差；但朱瑞华等［32］进一步加入气体热力学过程，并采用苏霍夫公式进行温度修正，计算结果与TGNET最接近。李欣泽［33］使用TGNET仿真软件计算出某管道的管存并确定了一些经典输量下管存的合理控制区间。丁媛等［34］以陕京天然气管道为例，使用 SPS（石油天然气长输管道模拟计算软件）进行瞬态模拟，由于瞬态模拟相较于稳态模拟考虑到了流量随时间的变化，减小了流动参数随时间变化的误差。使用仿真软件虽然能够较快地计算出管存，但若该仿真系统不是在线仿真系统，那么工况的变化将会对管存计算产生较大影响［33］，同时仿真软件需要人为设置管道的基础数据及边界条件，而往往设定值并不等于现场所采集的数据。

目前，国内外大型管道公司在管存运行管理方面也纷纷采用动态管存，如：国家管网集团北京油气调控中心将管道生产系统应用于实际运营中，通过 SCADA系统采集的实时数据计算管存以指导管道运行管理；哥伦比亚管道集团利用 IPS（基于工业互联网的智能化管道系统）将内外部数据充分整合，在此基础上搭建了工作管理系统、商业调度系统用以实时计算管存；TransCanada（加拿大横加管道公司）基于建设运营期的数据及跨业务信息的共享，实现管存计算及预测性分析功能。

由于管存计算的精确程度与所使用的计算方法密切相关，同时管道输气工况在气源供气量变化与用户需求波动的影响下常处于非稳定状态，而目前现场缺少国产自主研发的动态管存计算方法，因此，为适应天然气管道贸易精细化管理的要求，本文基于有限容积法开发一套具有自主知识产权的仿真算法，并以此为基础，进行天然气管道非稳定流动工况下的管存分析。在分析天然气管道稳态计算方法及常用的AGA8气体状态方程的基础上，推导了基于有限容积法离散后的天然气管道瞬态仿真模型，从理论上分析了稳态管存与动态管存计算方法的差异；使用现场实测数据验证了仿真模型的准确性，并对比分析了不同工况下管道稳态与动态管存差异，以指导管道运行方进行管存计算与输差分析。

1 稳态管存计算

1.1 稳态管存

管道某一时刻实际容纳的气体数量称为管存，通常用工程标准状况下的体积表示。在计算上有两种方法：稳态管存和动态管存。稳态管存即认为管道运行工况为稳定工况，根据QSY 05197—2019《油气管道输送损耗计算方法》中的天然气管道管存量计算方法，同时参考某国际天然气管道实际运行情况，可以用平均压力和平均温度表征管内气体当前状态，并以此为依据计算管存，具体计算方法如式（1）和式（2）。

式中：V—管道管存量，m3；p1—管道入口压力，Pa；p2—管道出口压力，Pa；pav—管道平均压力，Pa；V1—管道几何容积，m3；Z1—管道平均压力、平均温度下的气体压缩因子；T1—管道平均温度，K；Z0—工程标准状况下的气体压缩因子；T0—工程标准状况下的温度，取293.15 K；p0—工程标准状况下的压力，取101.3 kPa。

1.2 物性计算

天然气在不同压力及温度下的密度计算精度直接影响到天然气管道水力、热力计算及管存计算的准确性，因此确定合适的物性计算方法至关重要。参照GB/T 17747.2—2011《天然气压缩因子的计算第2部分：用摩尔组成进行计算》的推荐算法，同时，为了与某国际天然气管道（下文算例）物性计算方法保持一致，本文采用AGA8—92DC方程（简称AGA8方程）进行物性计算，在输入的参数误差满足标准且气体的热力状态在公式应用范围之内时，AGA8方程的不确定度小于 0.1%。AGA8方程于 1992年发表在 AGA（美国煤气协会）的 AGA8报告上，主要适用于绝对压力 0～12 MPa、温度263～338 K的管输天然气。该方程可输入多达 21种组分的摩尔比例进行计算，同时要求组分的摩尔分数大于0.000 05；对于H2S和CO2，适用的摩尔分数范围分别为0～0.02%，0～30%。

2 基于瞬态仿真的动态管存计算

2.1 动态管存

在天然气管道的运行中，常常伴随着运行工况的调整，如管道的进气量发生变化、压缩机的启停、调节阀的动作、沿线的下载气量发生变化等，这些都会导致管道处于非稳态，在这种情况下沿线的压力分布与稳态时的并不一致，因而不能采用平均压力公式（1）计算管道的压力。动态管存计算则是根据管道的基础数据，以管道两端的实测压力、温度和流量为边界条件，实时计算出沿线各个位置的压力、温度、流速和密度等参数，再根据管存计算公式（3），计算出管道的管存。管道各个位置的压力、温度、流速和密度等气体状态参数需要通过流动仿真获得。

式中：ρ—管段实际压力、温度下天然气的密度，kg/m3；0ρ—工程标准状况下天然气的密度，kg/m3；Vd—管段的设计管容量，m3；Vs—管段在标准参比条件下的管存量，m3。

2.2 基本模型

天然气管道流动仿真就是采用离散手段将无法直接得到解析解的描述天然气管道流动的偏微分方程组转化为一系列离散代数方程组，借助计算机完成数值求解，再通过对数值结果的分析再现管道内天然气的不稳定流动状态。用来描述天然气在管道内流动状态的控制方程称为管道的数学模型，是由连续性方程、动量方程和能量方程构成守恒型的天然气管流的基本控制方程。其中前两者可称为水力方程，后者称为热力方程。

式中：A—管道横截面积，m2；t—时间，s；w—气体流速，m/s；x—管道轴向长度，m；p—管道中气体的压力，Pa；D—管道内径，m；θ—管道倾角，°；λ—摩阻系数，无量纲；g—重力加速度，m/s2；e—气体的内能，J/kg；h—气体的焓值，J/（kg·K）；s—管道高程，m；Qq—天然气与环境之间的换热量，W/m3。

2.3 基于有限容积法的方程离散

天然气管道瞬态仿真的主流数值算法主要有特征线法、有限差分法、有限元法、状态空间模型法及有限容积法等。

特征线法其特点为显式求解，同时稳定性、适应性较高；但求解步长受条件限制，在快瞬变流研究中较为广泛［35］。现较常用的有限差分法为隐式有限差分法，该方法虽然求解速度很慢，但求解步长不受条件限制，同时不受物性及状态的影响，在长时间瞬变流的求解中应用较为广泛［36］。有限元法相较于有限差分法计算精度有所提高，同时求解速度也较快，在求解大型复杂管网时也具有一定优势［37］。状态空间模型法计算速度相较于有限差分法大幅提高，这是由于该方法基于传递函数直接表征了管道出入口的流动关系，但对于管道瞬时沿线压力流量难以求解［38］。有限容积法因其离散格式具有稳定、精准等特点，同时求解步长不受限制，因此在快、慢瞬变流动方面都有广泛的应用［39］。

基于以上对天然气瞬态仿真方法的调研，本文采用交错网格下的有限容积法进行流动仿真。有限容积法的基本思想是将问题的求解域离散化成一系列小区域，对每个小区域进行积分，得到关于未知函数的代数方程，联立求解这些代数方程就可得到未知函数在节点的参数值。

为了方便积分，需要对能量方程（6）进行转换。根据热力学定律有如下两式：

式中：u—气体黏度，Pa·s；T—管道内气体温度，K；Cp—气体的定压比热容J/（kg·K）。

将式（4）、（5）、（7）和（8）带入式（6），整理可得：

式（4）、式（5）和式（9）组成了天然气管道仿真的控制方程组。在交错网格下采用有限容积法对控制方程组进行离散（见图1）。主变量中速度变量位于网格的边缘，用虚直线标记；压力和温度变量位于网格的中心，用黑点标记。同时，密度、热容、声速等物性参数也位于网格的中心。

图1 管道网格划分示意图

对于计算所需，但该位置并未定义的变量，采用一阶迎风格式—公式（10）进行确定，如网格边缘的压力变量或者网格中心的速度变量。式（10）中，i表示第i个网格节点。

采用上述变量定义方法，动量方程式（5）可在动量控制体中离散成公式（11），式（11）中：上标“#”代表变量的当前值，“*”代表该变量下一次迭代中待求解的值。同理，连续性方程和能量方程可以在质量控制体中分别离散成公式（12）、公式（13）。

该节中，采用SIMPLE（压力耦合方程组的半隐式方法）求解天然气管道瞬态仿真中的水力、热力状态［40］。将下一时步速度、压力和密度的值写成当前值加修正量的形式，如式（14）～式（16）。

式中：w′，p′和ρ′分别为速度、压力和密度的修正量。将式（14）～式（16）带入式（11），并忽略二阶修正项，可得压力修正方程：

该方法的整体仿真流程如图2所示。仿真开始时，根据稳态计算结果对网格上各点参数赋值。接着根据离散后的动量方程式（11）求得一个新的速度值w。这个速度值不能匹配当前的压力值，因此根据式（17）求解压力修正量p′。继而根据式（16）求解密度修正量ρ′。至此，压力、速度和密度均进行了更新。这个循环修正的步骤迭代到压力修正量足够小才会结束。之后根据式（13）求解各网格点的温度值，这样这一时步的各点的压力、流量和温度均可获得，接着求解下一步，直到整个仿真计算完成。

图2 交错网格下的有限容积法计算过程

无论是稳态管存算法还是动态算法，其计算过程都是先设法获得管段中天然气的压力和温度，然后再根据该压力和温度进行其他计算，因此，确定压力和温度这一步骤十分关键。稳态算法利用天然气管道处于稳定流动时压力、温度的分布特点，给出了一种基于管道两端的压力、温度计算平均压力、温度的方法；动态管存算法则摒弃了这种思路，将管段划分为长度为1 km左右的小管段，再根据管道两段实测的压力、温度，仿真出每一小管段的压力、温度，最后再进行管存计算。显然后者管道划分更加精细，且对于非稳定工况下的管存描述更加准确。

3 应用实例

以某国际天然气管道 C1至 C2站间管道为例（管道参数如表1所示，气体组成如表2所示）。通过对比SCADA系统测量数据与仿真数据，验证基于有限容积法的天然气管网仿真方法的准确性。同时，针对此管道进行不同压力波动下的管存分析，用于分析管内流动的不稳定性对于管存计算的影响。

表1 管道参数

表2 气体组成

3.1 模型验证

通过SCADA系统获取该管道3 d的运行数据，包括C1站出站压力、流量与温度，以及C2站进站压力与温度。将C1站出站压力、C2站进站压力，C1站出站温度分别作为仿真的管道进口压力、出口压力及入口温度，仿真距离步长为1 km，时间步长为60 s，仿真160 h管内流动情况（见图3），共耗时3 min 44 s。将仿真所得C1站出站流量即管道入口流量与 C1站实测流量进行对比，平均相对误差为1.85%；同理，将仿真所得C2站进站温度即管道出口温度与 C2站实测温度进行对比，平均相对误差为0.92%（见图4）。以上结果验证了基于有限容积法的天然气管网仿真方法计算速度快且精度高，可以满足工业需求。

图3 仿真水力、热力边界

图4 仿真结果对比

3.2 动态管存分析

天然气管道内流动的不稳定性可以用压力变化程度进行表征。为了定量分析压力变化程度对于管存计算的影响，本文采用管道出口为定压，入口压力呈阶跃变化或周期变化两类算例进行仿真实验。对比入口压力波动，出口压力保持定值下，不同阶跃时长或不同周期下的动态管存与稳态管存差异。初始状态为入口压力8.5 MPa，入口温度45 ℃，出口压力6.5 MPa的稳定状态。

3.2.1 入口压力阶跃变化

压力阶跃变化即压力从一个稳定值线性增长到另一个稳定值，本算例中其数学模型为：

式中：Astep—阶跃振幅，MPa；Tstep—阶跃时长，s。Astep与Tstep共同决定了阶跃变化的形状，因此各设计3种阶跃幅度与时长，组合9种阶跃变化（见图5～图7）进行仿真实验，每次仿真时间为1 000 min，得到9种入口压力阶跃变化情况下的稳态与动态管存计算结果（见图8～图9）。

图5 入口压力阶跃变化下（阶跃时长6 h）稳态管存与动态管存对比

图6 入口压力阶跃变化下（阶跃时长2 h）稳态管存与动态管存对比

图7 入口压力阶跃变化下（阶跃时长0.5 h）稳态管存与动态管存对比

图8 入口压力阶跃变化下管存偏差对比

图9 入口压力阶跃变化下管存相对管存差对比

由图5可以看出，两类管存均是先上升，后平稳，总体趋势一致，但是稳态管存上升速度比动态管存快，且流动稳定后，稳态管存均较动态管存低。同时，稳态管存变化过程在t=Tstep时刻存在一个明显的拐点，该拐点之前由于入口压力快速增长，因此平均压力也快速增长，稳态管存也变化迅速；拐点之后，进出口压力均不变，但温度场并未达到稳定，出口温度缓慢下降，造成平均温度减小，稳态管存缓慢增长直至稳定。由此可知，由于稳态管存计算方法仅仅与进出口两个点压力与温度有关，因此阶跃变化造成的入口压力变化会立刻影响到平均压力的计算，而实际情况是入口压力的上升是管道充装气体造成的，而数百千米的管道充装是需要时间的，因此实际管存变化应该如动态管存所示，平滑变化，且当入口压力阶跃结束后，仍维持一段时间的管存上升状态。此外，流动稳定后，高压状态下的稳态管存与动态管存偏差变小，说明高压稳定状态下的平均法计算管存是相对有效的；相同阶跃振幅下，阶跃时长越短，管存上升阶段（0～Tstep）稳态管存变化斜率较动态管存变化斜率更大，说明在管内不稳定流动刚发生时段进行管存稳态计算会造成较大的误差，且振幅越大，误差越大。

在输差分析过程中，管存所计算的往往是相对量，即，终了状态管内管存量与初始时刻之差。因此本文对比了稳态与动态管存差之间的偏差，即：

式中：ΔLP—稳态管存差与动态管存差之间的偏差（简称管存偏差），104m3；dLsP—当前时刻稳态管存与初始时刻偏差，104m3；dLPd—当前时刻动态管存与初始时刻偏差，104m3。

由图6（a）可以看出，管存偏差变化趋势均为先增加后减小再平稳，且增加减小的速率均是逐渐变缓，拐点出现在t=Tstep时刻。分析原因可知，0～Tstep时间段内，上升速率大于动态管存，但稳态管存接近线性上升，而动态管存呈抛物线增长，因此管存偏差呈抛物线增长；当t＞Tstep，稳态管存上升速度迅速下降，低于动态管存，因此管存偏差又下降，直到流动稳定后平稳。相同阶跃时长下，振幅越大，管存偏差越大；相同阶跃振幅下，时长越长，最大管存偏差越小，但流动稳定后管存偏差相同。

以上为定性分析，以下从定量的角度分析管存计算方法对于输差计算的影响程度。参考相对输差计算公式［14］，引入相对管存差参数RLP，RLP越大，由管存计算不准导致的输差越大，反之亦然。

式中：RLP—相对管存差，无量纲；Qall—0～t时间段内的总进气量，104m3；L—输气管道总长度，km；l—某段管道长度，km。算例中l=202.8 km，L=524.4 km，因此RLP=2.59 △RLP/Qall，计算得到9种阶跃变化下的RLP（见图6（b））。从理论上分析，随时间跨度增大，输气量增大，由于管存计算偏差带来的输差会逐渐减小，这与图6（b）相符合。除最初几分钟外，RLP随着时间增加而降低。本文以3%为阈值，当RLP＜3%时，认为管存偏差对输差影响可以忽略不计，统计阈值时间点Ts，即t＞Ts，RLP＜3%，可以得到入口压力阶跃变化下阈值时间点（见表3）。此外，将阈值时间点与阶跃时间相减可得额外稳定时间，用于描述阶跃变化停止后到管存偏差对输出影响可以忽略不计的时间。显然，阶跃时长越短波动越快，所需额外稳定时间也越长。如，在阶跃振幅都为0.5 MPa时，阶跃时长为6 h时额外稳定时间最短。阶跃时长为 0.5 h时稳定时间最长同时，阶跃振幅越大，所需额外稳定时间也越长。如，在阶跃时长都为6 h时，额外稳定时间随着振幅的增加而延长。

表3 入口压力阶跃变化下阈值时间点

3.2.2 入口压力周期变化

实际天然气管道运行工况受上下游供需影响，而供需受温度等周期变化因素影响而存在一定周期性。因此本文设计正弦周期变化的算例，用于分析不同周期与振幅的正弦入口压力变化下，稳态管存与动态管存计算结果偏差，其数学模型如下式：

式中：Asin—正弦振幅，MPa；Tsin—正弦周期，s。Asin与Tsin共同决定了正弦变化的形状，因此各设计3种正弦振幅与周期，共组合9种正弦变化（见图10～图12），进行仿真实验。

图10 入口压力周期变化下（正弦周期24 h）稳态管存与动态管存对比

图11 入口压力周期变化下（正弦周期12 h）稳态管存与动态管存对比

图12 入口压力周期变化下（正弦周期2 h）稳态管存与动态管存对比

总体来看，两类管存变化趋势与入口压力变化相似，均呈现正弦变化，但是动态管存计算结果相较于稳态管存存在一定程度的相位差与振幅差。相同Asin、不同Tsin下，稳态管存的正弦变化振幅基本不变，周期也与Tsin保持一致，这是不合理的。在入口压力持续变化情况下，管内始终处于不稳定流动状态，入口压力波动传递到200 km管长的各处是需要时间的，造成的管存变化必然滞后。动态管存计算方法则能较好体现这一特性，管存的变化较之入口压力变化存在一定程度的时间滞后。此外，随着Tsin变小，入口压力波动变大，管内气体无法充装到稳态管存量的峰值就下降，因此动态管存的峰值随着Tsin变小而变小，这与现场实际情况相符。

类似于阶跃变化算例，周期变化算例也对稳态与动态管存差之间的偏差进行了分析（见图 13）。入口压力变化振幅越大、周期越短，管存偏差越大，也就是管内流动波动越大，稳态管存计算方式越不适用，这与直观理解是相符合的。此外，关于RLP的结果（见图14、表4）所示，其结果与阶跃变化下的结果略有不同。周期变化下的RLP随时间围绕 0振荡衰减，且短周期高振幅的周期波动所需要的阈值时间点很长，如，周期2 h、振幅1.5 MPa算例中，阈值时间点为4 232 min，需要接近3 d的时间才可以忽略管存计算对于输差计算的影响。

表4 入口压力周期变化下阈值时间点

图13 入口压力周期变化下管存偏差对比

图14 入口压力周期变化下相对管存差对比

3.3 小结

通过与现场实测数据对比，验证了基于有限容积法在天然气管网仿真算法方面的可靠性，同时通过对于入口压力阶跃变化与周期变化两种情况的稳态管存与动态管存计算方法对比可知，动态管存计算方法的适用性与精确性远高于稳态管存计算方法。对于管道运行方来说，应改用动态管存计算方式进行管存管理与输差分析，以减少由于管存偏差而导致的经济损失。同时以下3点还需进一步优化。

（1）现场数据存在大量噪声与异常，若直接进行计算会产生计算偏差，首先应进行数据清洗，通过提高数据可靠性有效降低计算偏差。

（2）需进一步研究SCADA系统和计算程序之间的通信问题，这样在保证计算精确度的同时也提高了运行管理的效率。

（3）计算模型的精确程度直接影响着现场运行管理，因此还需要基于运行数据不断地修正计算模型，以提高模型的自适应性来感知现场工况的变化。

4 结论

精细化管理是天然气管道运行的重要要求，管存计算也逐渐成为各管道公司日常运行管理内容。管存计算的精确程度不仅涉及输差问题，同时也会对于管道日常运行、储气调峰及应急供气产生影响。本文从天然气管道瞬态流动仿真出发，提出一套完整的动态管存计算方法。该方法以有限容积法为核心，求解流动三大基本方程，并结合AGA8气体状态方程，精确计算天然气管道非稳定流动状态下沿线压力温度密度分布，进而获得管内管存精确变化情况。本文用现场数据对该仿真方法进行了验证，并对比了现行稳态管存与本文提出的动态管存计算方法在不同阶跃、周期变化下的偏差，体现了动态管存计算方法的优势。此外，随着天然气管道运行数字化、智能化，可以预见将动态管存计算方法与在线仿真相结合，可以获得实时管内管存变化情况，为运行方提供更精确的决策支持信息。