天然气管网仿真分层计算方法

阎 涛，刘天尧，王雪莉，薛向东，2，康 阳，朱 峰

1.国家管网集团科学技术研究总院分公司，河北廊坊 065000

2.西安交通大学软件学院，陕西西安 710049

随着经济的持续发展，我国对能源的需求越来越大[1]。由于天然气具有高热值、低碳、清洁等优点[2-3]，其消费量逐年增长。预计到2030 年，全国天然气消费量将达到（5 500～6 000）×108m3[4]。出于经济原因，天然气通常通过管道进行大规模输送[5-6]。通过对天然气管网进行仿真，可以明确天然气在管网内部流动过程中的水热力变化情况，为天然气管网的优化、管理、调度及运行提供技术支持[7-10]。截至2021年，我国主干天然气管道总里程达到11.6×104km[11]。天然气管网规模的不断扩大对管网仿真提出了越来越高的要求。对于大规模天然气管网，采用隐式差分法离散水热力系统得到的代数方程组规模较大。受限于代数方程组巨大的规模，LU 分解等直接解法的[12]求解速度慢。童睿康等[13]采用迭代法替代高斯消元法来求解水热力方程组，有效提高了求解速度。然而，使用迭代法仍要对整个方程组进行求解操作，这对大规模管网来说仍是个较大的挑战。为提高天然气管网的求解速度与效率，提出了一种管网分层求解的算法。该算法将管网划分为不同层次的子管网，自下而上传递变量关系式，自上而下传递变量值，最终实现整个管网仿真过程的快速求解。

1 管网数学模型

天然气管网由管道、阀门、压缩机等元件组成，因此描述气体流经这些元件行为的数学模型就组成天然气管网的数学模型。

1.1 管道元件数学模型

天然气管道的数学模型由连续性方程、动量方程和能量方程共同组成。其中，连续性方程和动量方程被称为管道的水力方程，能量方程则被称为管道的热力方程[14]。采用质量流量m、压力p、温度T作为描述管道内部水力和热力过程基础变量[15]的方程形式如下式所示：

连续性方程：

动量方程：

能量方程：

式中：p为管道内气体压力，Pa；t为时间变量，s；T为管道内气体的温度，K；A为管道横截面积，m2；ρ为天然气密度，kg/m3；m为管道内气体的质量流量，kg/s；z为轴向空间变量，m；λ为水力摩阻系数；g为重力加速度，m/s2；θ为管道倾斜角，rad；cv为气体的定容比热容，J/（kg·K）；u为管道内气体流速，m/s；D为管道内径，m；Kt为管道内气体与环境之间的总传热系数，W/（m2·K）；Tg为管道周围土壤的温度，K。

其中，质量流量m和压力p称为水力变量，温度T称为热力变量。

为了加快管网仿真过程求解速度，可以对上述方程进行线性化处理。线性化的统一表达式[16]如下：

式中各参数的具体形式如表1所示。

表1 管道数学模型参数

采用Implicit Cell Centered Method[17]方法对线性化后的水力方程进行离散，整理后可以得到如下方程：

线性化热力方程的时间项采用向前差分格式，对流项采用一阶迎风格式进行离散，离散结果整理后如下：

式中：UPi、CEi、DWi为根据上一时层计算得到的热力系数。

1.2 非管道元件数学模型

天然气管网元件可以分为管道元件和非管道元件，其中非管道元件包括压缩机、阀门、气源、节点等。非管道元件具体的数学模型可参考相关文献[18]。

2 管网求解策略

为加快天然气水热力仿真的速度，采用水热力解耦[19]的求解策略。如图1所示，该方法将天然气管网的水热力过程解耦为水力过程和热力过程再分别求解，能够在保证仿真精度的前提下提高仿真速度。

图1 天然气管网水热力解耦算法

2.1 分层算法计算简介

如图2所示，直接求解算法是直接采用LU分解算法对整个天然气管网的水力（热力）方程组进行求解。由于天然气管网元件众多且管道离散节点较多，该方程组规模较大，直接求解算法较为耗时。

图2 分层算法原理

针对这一问题提出了一种分层算法。该算法主要针对不同的子管网集合进行计算，因此在应用分层算法之前需要将管网按照一定的标准分解为若干子管网的集合。对于常见的燃气管网或者长输天然气管网，可以按照调压站（调压阀）或增压站（压缩机）所在的位置对管网进行分解操作，从而可以形成不同层级的子管网。需要指出的是，本文提出的分层算法包括图2 的步骤2 和步骤3，但是步骤1并不属于分层算法的内容。

将管网分为若干个子管网之后，将其组织为层级结构，之后便可以采用分层算法，开始通过LU 分解方法对各个子管网进行单独求解，求解过程中采用CPU 并行技术[20-21]将同一层级的子管网同时进行计算，从而加快仿真求解速度。

2.2 分层算法求解水热力过程

图3 代表由a层b个子管网通过拓扑连接形成的天然气管网。子管网内部元件之间的拓扑结构并不会影响分层算法的应用，因此图3并未显示各个子管网内部的拓扑结构。管网分层操作使得同一个物理节点同时属于两个子管网，通过虚线来表示这种子管网间的拓扑连接。

图3 具有多层结构的天然气管网

如图4 所示，某层的子管网α 与m个子管网连接，其中有n个顶层子管网和m-n个底层子管网。需要说明的是，顶层子管网和底层子管网都是相对本层子管网而言的。

图4 与外界连接的子管网

假设图4中的子管网α内部有Np个管道、Nv个阀门、Nc个压缩机、Ne个气源以及Nd个节点。在水热力求解过程中，为了保证管道内部变量的求解精度，通常将管道离散为较多的管段。此处为了简化说明过程，子管网内部的管道没有进行离散操作，但这并不会影响分层算法的具体实施流程。

管道、阀门和压缩机都具有起点和终点两个端点，气源只有一个端点，而每个端点上都具有压力和流量两个水力变量，因此子管网α内部管道元件、阀门元件、压缩机元件、气源元件端点处的水力变量数目共计4Np+4Nv+4Nc+2Ne。

管道、阀门和压缩机都具有2 个特性方程，而气源具有1 个特性方程为已知的边界条件值，因此，子管网内部管道水力方程，阀门、压缩机和气源的水力特性方程共计2Np+2Nv+2Nc+Ne个。根据图论知识可知，与m个外界子管网连接的子管网内部节点上流量平衡方程和压力平衡方程数目为2Np+2Nv+2Nc+Ne-m。至此，方程数目为4Np+4Nv+4Nc+2Ne-m。

子管网α内部水力变量数目较水力方程数目多m，方程不封闭。为了使水力方程组达到封闭的条件，需要补充m个方程。

如图4所示，本层子管网与m-n个底层子管网连接。这些底层子管网向子管网α传入m-n个补充方程，子管网中水力方程组数目变为4Np+4Nv+4Nc+2Ne-n。该补充方程为连接点上两个水力变量之间的关系式，即连接点上质量流量m和压力p之间的关系式。

子管网α 通过n个连接点与顶层子管网连接，选用该n个连接点的质量流量m（或压力p）作为自由变量v，可将管网内部所有水力变量表示为关于该n个自由变量（v1，v2，…，vn）的关系式。

令[v1，v2，…，vn]=[1，0，…，0，…，0 ]，代表关于自由变量（v1，v2，…，vn）的n个方程，其中v1= 1，v2～n= 0。

将该n个方程作为4Np+4Nv+4Nc+2Ne-n个方程的补充方程，通过求解这些方程组可得解向量D1。求解方程组可以采用LU 分解等常见的方程组求解算法。

将关于自由变量（v1，v2，…，vn）的n个方程中的各个变量依次等于1，重复上述过程，于是得到一系列的解向量D2，D3，Dn。

最后令[v1,v2,…,vk]=[0,0,…,0,…,0 ]，求解得到Dn+1。

于是可得子管网α内部所有水力变量关于选定的n个自由变量（v1，v2，…，vn）表达式：

根据上式可得连接点上选定的自由变量v与非自由变量之间关系式：

上式中c表示与子管网α 连接的子管网编号。若自由变量v代表压力p，非自由变量vˉ则代表流量m；若自由变量v代表流量m，非自由变量vˉ则代表压力p。

将上述n个关系式作为补充方程传入顶层子管网中，正如之前由底层子管网传上来的m-n个补充方程。

如图5所示，底层的子管网不断向顶层的子管网传递补充方程，直至最顶层的子管网。最顶层的子管网待底层子管网传入补充方程后，水力方程组闭合。求解该方程组便可得到最顶层的子管网内部所有变量的值，包括最顶层的子管网与其底层子管网连接点上自由变量的值。值得一提的是，分层算法使得每层之间的各个子网的求解互相解耦，因此可利用并行计算等方式实现计算速度的提升。

图5 分层算法示意

将这些自由变量的值传入底层子管网，根据前述子管网内部计算得到关系式，得到子管网内部所有变量的值。不断重复该过程至最底层管网，从而得到整个管网水力变量值，完成管网水力过程的求解。下面以一个具有详细拓扑结构的天然气管网为例展示分层算法求解水力仿真的过程。

如图6 所示的天然气管网具有9 个气源，S1 为进气气源（源点），S2～S9 为出气气源（汇点），其中进气气源可以代表储气库或上游管网来气。该天然气管网通过阀门V1～V6 调节压力。实际管网中存在多种环状管网，不能一一列举，因此图中只采用三角形环状管网，但这并不意味着分层算法只能应用于三角形环状结构的管网。

图6 具有多个调压阀的天然气管网

按照调压站（阀）的位置将管网分为如图7所示的具有3 个不同层级的子管网，其中最底层有4个子管网（3-1，3-2，3-3，3-4），中间层有2 个子管网（2-1，2-2），最顶层有1个子管网（1-1）。图中红色虚线代表分层之后不同子管网之间的拓扑连接关系。由于管网在阀门处进行分层，因此选用分界处的阀门变量作为自由变量。

图7 分层后的天然气管网拓扑图

子管网3-1 共有4 条管道，1 个阀门和2 个汇点。依旧假设管道没有离散，只有起点和终点两个节点。子管网3-1内部的水力变量和方程数目如表2所示。

表2 子管网3-1水力变量与方程数目

子管网3-1 与外界有1 个连接点，从表2 可知其水力变量数目较方程数目多1。为了使方程组闭合，需补充1 个关于自由变量p5的方程p5=1，并且其余方程右边的常数项设置为0。采用LU 分解算法求解闭合的方程组可得到解向量A1。同样，在原先缺1 个方程的方程组基础上再添加p5=0 的方程，其余方程右边的常数项也设置为0，采用LU分解算法求解闭合的方程组可以得到解向量B1。因此可以得到如下关系式：

其中，A51和B51代表解向量A1和B1中变量m5对应的系数。

子管网3-1 和3-2 将分别将表达式m5=p5A51+B51和m9=p9A92+B92传递给顶层子管网2-1。子管网2-1 与外界有3 个连接点，因此其需要3 个方程才能达到方程组闭合的条件，底层子管网已经传递了2 个方程，因此其只需要补充1 个方程即可。接着重复子管网3-1的求解流程即可。

表3 是不同层的子管网应用分层算法的过程中各个连接点选取的自由变量以及产生的一些方程。

表3 不同层的子管网应用分层算法时自由变量与方程

采用分层算法对天然气管网进行热力求解与上述水力求解类似，此处不再赘述。

3 案例分析

某天然气管网具有如图8 所示的复杂拓扑结构，该管网具有105 条管道（见表4）、3 台压缩机、7个阀门以及32个气源。

图8 具有复杂拓扑结构的某天然气管网

表4 管道长度

所有管道的规格均为D720 mm×10 mm，出气气源（S1～S14，S16～S32）都采用流量边界（见表5），进气气源（S15）采用压力边界（3 MPa），压缩机（C1～C3）具有相同的压比-流量曲线（见表6）。

表5 各个出气气源（汇点）的流量

表6 压缩机的流量压比数据

管网t=0时刻已经处于稳态。当t=2 h时，阀门V1关闭，管网停止向出气气源S2供气。

按照管网中调压阀V2、V3、V6、V5和增压站C2、C3 的位置将管网分为如图9 所示的不同层级子管网集合。

图9 分区后的天然气管网拓扑结构

3.1 准确性对比

采用分层算法对上述管网算例进行求解。由于停止向气源S2 供气之后，管道P4、P5、P6 及气源S1、S3、S4 这些元件离关闭的阀门V1 较近，受到的影响最为明显，因此选择将管网停止向气源S2 供气前2 h 和后24 h 总共26 h 的这些元件的压力或流量进行对比。同时为了保证对比的充分性，还将对比管网再次达到稳态后采用分层算法与不采用分层算法（即直接求解算法）计算的不同位置处管道（P7、P8、P9、P10、P11、P12）、气源（S5、S9、S19、S21、S26）、压缩机（C1、C3、C4）的流量、压力以及温度。

需要说明的是，两种算法对管网仿真过程中出现的一些小方程组都采用LU 分解法进行求解。直接求解算法与分层算法的水力和热力仿真时间步长为20 s，水力和热力计算的空间步长都为1 km。

如图10～图14 所示，无论是关阀事件发生之后一段时间内阀门附近受到影响最大元件的流量、压力和温度对比情况，还是关闭阀门之后天然气管网再次达到稳态时不同管道、气源和压缩机的流量、压力和温度对比情况都表明，分层算法计算得到计算结果与直接求解的计算结果完全一致，这说明分层算法具有和直接求解算法同等的计算精度，分层算法的正确性和准确性得到验证。

图10 管道P4、P5、P6起点流量随时间变化情况

图11 管道P4、P5、P6起点压力随时间变化情况

图12 管道P4、P5、P6起点温度随时间变化情况

图13 气源S1、S3、S4压力随时间变化情况

图14 气源S1、S3、S4温度随时间变化情况

3.2 求解速度对比

对同在一层的各个子网，可以利用CPU 并行计算实现子管网计算过程的加速。本实验基于AMD Ryzen 72700x 的CPU 平台进行测试，表7 是分层算法与直接求解算法在不同管道总长度的管网算例上仿真耗时的对比情况。

表7 两种算法仿真耗时对比

由表7 可知，在不同规模管网的仿真过程中，分层算法都具有较快的仿真速度，加速比平均在2倍左右。值得一提的是，本案例中仅对比了同样拓扑结构不同管网规模下的加速效果，可知对其他拓扑结构的管网，当每层管网的独立子网数量越多时加速效果会更加明显。

4 结论

针对大规模天然气管网水热力模型计算时直接求解算法求解速度慢的问题，提出了一种将管网分层的算法。该算法将管网分为若干层的子管网，通过自下而上的补充方程传递和自上而下的值传递过程完成管网方程组的求解。通过天然气管网向汇点供气的算例来说明分层算法应用于管网水热力仿真过程的计算精度和速度。算例的计算结果表明：本算法的应用不会改变模型求解结果的精度，能够确保基础模型的准确性；本算法可以结合并行计算实现计算加速，对本案例中不同规模的天然气管网仿真计算，平均可达2倍左右的加速比。