0 引言

随着天然气这一清洁能源消耗量的增加，天然气管网的规模在不断扩大，对输气管道内稳态流动准确模拟的需求愈加迫切，这对天然气输气管道设计及管理具有重要意义^[1]。

输气管道的稳态仿真就是对定常流动的连续方程、动量方程、能量方程和状态方程等进行求解的过程^[2-3]。等温输气管道水力计算公式是基于等温输气管道模型，通过对连续方程和动量方程进行积分得到的，没有考虑管道温度的变化，如果管道温度变化较大，仿真计算会存在较大误差。笔者在此基础上引入能量方程，并提出分段积分法，使其可以对非等温输气管道进行稳态仿真。以川气东送干线输气管道为例进行仿真计算，并与TGNET软件和龙格库塔方法仿真结果进行对比。此外，为了方便管道运行能耗最优化方案计算，对已知管道末端压力及温度条件下的稳态仿真方法进行了阐述，并给出了相应的计算实例。

1 管道稳态流动控制方程

气体在输气管道中的流动可以看作一元流动，稳态流动与时间无关，略去控制方程中的时间相关项，沿流动方向输气管道稳态流动模型如下。

连续方程：

(1)

动量方程：

(2)

能量方程：

(3)

状态方程：

(4)

焓方程：

(5)

式中：ρ为气体密度，kg/m³；w为气体流速，m/s；p为气体压力，Pa；g为重力加速度，m/s²；θ为管道与水平面之间的夹角，rad；λ为管道的摩阻系数；d为管道内径，m；Q为单位质量气体的吸热或放热量，J/kg，；h为气体的焓，J/kg；s为管道末端与始端高程差，m；T为气体温度，K。

在已知混合气体组分和混合气体的压力、温度条件下，可以由状态方程求解气体的密度、压缩因子和比热容等，例如SRK状态方程、PR状态方程和BWRS状态方程等^[4-6]。笔者在研究中采用BWRS方程计算，其中管道摩阻系数λ采用科尔布鲁克公式计算^[7]。

2 方程求解方法 2.1 龙格库塔法

以管道沿线的温度、压力和流量为未知量，可以将式(1)~式(5)整理成以下的形式：

(6)

其中：a₁₁=0, a₁₂=0, a₁₃=1, b₁=0，a₂₁=

式中：m为管道内气体质量流量，kg/s；A为管道横截面积，m²；T_e为管道埋深处土壤温度，K。

式(6)是一组常微分方程组，管道起始端的压力p₀、温度T₀和流量m₀是已知条件。文献[8-11]采用了四阶龙格库塔法对微分方程组进行求解，需要将输气管道划分成N段微元管道。这里以第i段微元管道计算为例进行说明。

(1) 利用第i段管道的初始端压力p_i、温度T_i和流量m_i，可以求解出方程组(6)中的系数，进而求解出(dp/dx)₁、(dT/dx)₁和(dm/dx)₁，然后按照下式得到p_i¹、T_i¹和m_i¹(式中Δx为每段微元管道的长度)。

(7)

(2) 利用新得到的p_i¹、T_i¹和m_i¹求解出(dp/dx)₂、(dT/dx)₂和(dm/dx)₂，然后计算p_i²、T_i²和m_i²。

(8)

(3) 利用新得到的p_i²、T_i²和m_i²再求解出(dp/dx)₃、(dT/dx)₃和(dm/dx)₃，然后计算p_i³、T_i³和m_i³。

(9)

(4) 利用新得到的p_i³、T_i³和m_i³求解得到(dp/dx)₄、(dT/dx)₄和(dm/dx)₄。然后计算该段管道末端的压力p_i+1、温度T_i+1和流量m_i+1。龙格库塔方法的计算公式如下：

(10)

从管道起始端开始，反复使用公式(10)，即可完成从管道起始端到末端的压力、温度和流量计算。这种求解常微分方程组的方法对非等温输气管道的仿真也有很高的准确性。

2.2 分段积分法

首先将非等温输气管道划分成N段管道，当N足够大，每段管道内气体流动过程的温度变化则很小，可以用管道的平均温度作为管道水力计算的参考。首先计算每段管道的平均温度，然后采用等温输气管道水力计算公式求得管道的压力分布，按照气体流动方向，依次对N段管道进行计算，一直求解到管道末端。

管道起始端的压力p₀、温度T₀和流量m₀是已知条件，这里以第1段为例介绍管道的水力和热力计算方程。对管道连续方程、动量方程和状态方程联立积分得到等温输气管道水力计算公式^[12-13]，已知管道起始端压力和气体流量就可以求得末端压力：

(11)

对管道能量方程和焓方程进行积分得到输气管道热力计算公式，已知管道起始端气体温度和沿线地温，就可以求得管道内气体温度(当x=L时，得到管道末段温度)。

(12)

按照每段管道的全长，对压力积分，求得每段管道平均压力：

(13)

若输气管道温降曲线与沿线坐标所包围的面积，同某一恒定温度与沿线坐标所包围的面积相等时，称这一恒定温度为管道的平均温度，代入式(12)得到平均温度表达式：

(14)

式中：q为气体体积流量，m³/s；C为计算常数；a=2g/ZRT；Z为气体压缩因子；R为气体常数，m²/(s²·K)；Δ为气体的相对密度(标况下气体与干空气的密度之比)；L为管道长度，m；s_i为相对于第i段管道起点的高程，m；L_i为第i管道的长度，m；p₀、p₁分别为管道起始端和末端压力，Pa；T为管道长度x处气体温度，K；T₀为管道起始端气体温度，K；α=Kπd/mc_p；K为管道总传热系数，W/(m²·K)；c_p为气体质量定压热容，J/(kg·K)；D_i为焦耳-汤姆逊系数, K/Pa；p_cp为管道平均压力，Pa；T_cp为管道平均温度，K。

式(11)~式(14)中的部分系数采用状态方程计算得到。对第i段管道仿真计算的具体步骤：

(1) 已知第i段管道起始端的压力p_i、温度T_i和流量m_i。假定第i段管道末端压力值p_i+1′，由式(13)和式(14)分别计算得到管道平均压力p_cp和平均温度T_cp；

(2) p_cp、T_cp分别作为系统压力和温度，计算得到气体的压缩因子和密度等参数，由式(11)得到新的管道末端压力值p_i+1″；

(3) 更新管道末端压力值，重复(1)、(2)步，直至前后2次计算得到的管道末端压力值之差满足精度要求，停止计算，得到管道末端压力p_i+1，由式(12)计算末端温度值T_i+1；

(4) 以p_i+1和T_i+1作为下段管道起始端压力和温度，对下段管道进行计算。

从管道起始端开始，对N段管道依次计算，就可以完成对管道的稳态仿真计算。

2.3 已知管道末端条件的仿真方法

管道的仿真计算给定的已知条件一般都是管道起始端的压力、温度和流量条件，由管道起始端向管道末端进行计算。而在管道运行方案优化计算中，常常遇到给定管道末端压力、温度和流量条件，需要求解管道起始端的压力、温度和流量情况。现对这种情况的计算过程进行说明。

针对龙格库塔方法，将式(7)~(10)做如下修改，即有式(15)~(18)。式中下标N代表管道末端条件，下标M代表管道起始端条件。

(15)

(16)

(17)

(18)

按照2.1节所述龙格库塔方法进行求解。

在使用分段积分法时，也需要对式(11)~式(14)进行相应修改，即有式(19)~式(22)。计算步骤与2.2节类似。

(19)

(20)

(21)

(22)

3 计算实例

以川气东送干线输气管道作为计算对象。川气东送管道西起川东北普光首站，东至上海末站，途经8省53个县市。管道全长2 229 km，包括1条干线、5条支线和1条专线。干线设计压力10 MPa，设计输量120×10⁸ m³/a。气源为普光净化厂来气、元坝气田气和涪陵页岩气混合气输送^[14]。管道干线包括普光首站、利川站和潜江站3个增压站，节点包括普光出、梁平、黄金、王场、利川和恩施等28个节点，管道包含普光出-梁平、梁平-黄金、黄金-王场、王场-利川入、利川出-恩施等27条输气管道。采用真实的气源信息和管道基础数据进行计算与分析，气源进气量、节点下气量和管道沿线地温采用某天的运行参数，压气站出口压力设置为8.000、8.300和7.500 MPa，出口温度为45.000 ℃。

TGNET是美国SSI公司的Pipeline Studio管道软件分析包中对输气管道进行工艺分析的软件，能够对输气管道中的单相气体进行稳态和动态模拟^[15]，是管道模拟方面的权威软件之一。笔者分别采用TGNET软件、龙格库塔法和分段积分法对川气管道进行仿真计算。为了方便对比，计算中取管道分段步长为相同值，这里把龙格库塔法和分段积分法的管道分段步长设为50 m。表 1中按照气体流动方向，列出了各个节点计算得到的压力和温度信息。从表中可以看出，笔者所用的2种方法进行管道稳态仿真的计算结果与商业软件TGNET的计算结果非常接近，细微的差别可能是TGNET软件与编程计算的气体参数不同带来的。分段积分法的计算结果更接近TGNET的结果，管道末端节点(上海站)的压力值误差最大，龙格库塔法的计算误差为1.847%，分段积分法的计算误差为1.681%。笔者认为分段积分法误差更小的原因在于管道的计算过程中对压力和温度不断进行迭代计算，修正了输气管道中流速和温度变化对于管道仿真计算的影响。

由于分段积分法和龙格库塔法是间断性的求解，会产生积累误差，所以输气管道的分段步长要尽量小，以减小积累误差。表 2中列出了管道分段步长从500~50 m这7种情况下节点压力值的最大误差，并比较了2种方法所需要的计算时间。从表中数据可知，当管道分段步长较大时，分段积分法的最大误差更大；当管道分段步长减小后，分段积分法的最大误差减小很快；当分段步长小于100 m后，分段积分法的最大误差基本稳定在1.670%左右，小于龙格库塔方法的计算误差1.847%。计算时间方面，相同分段步长条件下，分段积分法的计算时间小于龙格库塔方法，且分段积分方法的计算效率更高。

以利川出-恩施管道为例，在已知管道末端流量、压力及温度的条件下，进行管道仿真计算。气体由利川出站向恩施站流动，所以利川出站为管道起始端，恩施站为管道末端。利川出-恩施段标况流量2 275.834×10⁴ m³/d，地温18.800 ℃。采用恩施站计算结果作为已知条件，龙格库塔法已知恩施站压力7.991 MPa，温度29.387 ℃，采用式(15)和式(18)进行计算，求利川出站的压力和温度。分段积分法已知恩施站压力7.997 MPa，温度28.743 ℃，采用式(19)~式(22)进行计算，求利川出站的压力和温度。取管道的五等分点为观察节点，各节点求得的温度和压力结果见表 3。理论上正确的计算结果是利川出站压力8.300 MPa，温度45.000 ℃。由表 3可知2种方法的压力和温度沿管道变化趋势正确，利川出站的压力和温度计算正确，验证了仿真计算方法的正确性。

4 结论

提出管道稳态仿真的分段积分法，在等温输气管道水力计算公式中引入温度的计算，用来对非等温输气管道进行稳态仿真计算。以川气东送干线输气管道为求解对象，采用龙格库塔法和分段积分法进行计算，并对比TGNET的计算结果，得到以下结论：

(1) 龙格库塔法及分段积分法都有较高的计算准确性，可用于输气管道的稳态仿真计算。

(2) 当管道分段步长较大时，龙格库塔方法精度更高；管道分段步长较小时，分段积分法的精度更高。管道分段步长相同的条件下，分段积分法的计算效率更高。

(3) 这2种方法同样适用于已知管道末端压力、温度及流量条件下的稳态仿真计算。