企业官网建设流程全解析

简介：6节点天然气潮流计算程序是一套基于MATLAB的教学实例，面向能源、化工等专业初学天然气网络分析的学习者，用来计算6节点模型中压力、流量、存储量等关键参数，理解网络潮流分布与迭代求解逻辑。压缩包共2个文件，均为m格式，包含主程序与核心函数脚本，整体仅1KB，结构轻量、注释清晰，适合逐行研读和二次修改。目前已有499人学习。程序覆盖了天然气状态方程、节点能量守恒、管道阻力与压降计算，以及牛顿法等非线性方程组的迭代求解思路，同时展示从数据输入、核心计算到结果可视化的完整代码框架。借助此程序，学习者能更直观地将理论转化为实际可运行代码，也为后续扩展至更复杂的天然气管网分析打下基础。 很多人第一次听到“天然气潮流计算”这六个字，容易觉得是电力系统的人跑来抢饭碗。其实不管是电网还是天然气管网，只要管道一多、负荷一杂，靠手算或者凭经验拍脑袋都是要翻车的。我最近整理了一个6节点天然气潮流计算的MATLAB程序，把天然气门站到几个用户节点的稳态压力、流量分布一次性算清楚。这篇文章就把这个程序的建模思路、求解器设计、踩坑经验和结果验证完整拆开讲，适合正在做管网仿真、综合能源系统课设，或者刚接触气体网络计算的朋友直接照着复现。

这里说的“潮流计算”，本质上就是求解一组非线性方程：给定气源压力和节点用气负荷，反推出整个管网里每个节点的压力和每条管道的流量。6节点规模不大，但足够把环网、多负荷、压力降、回流这些核心现象都展示出来，比动辄几十上百节点的工业算例更容易看清楚算法本身。

1. 为什么6节点天然气潮流计算值得自己写一遍

接触过电力系统的人都知道，电网潮流是电力专业的基础课，有牛拉法、PQ分解法一大堆成熟工具。天然气系统这几年在综合能源、双碳背景下被反复提起，可很多人对气网的认知还停留在“一根管子送气”的阶段，遇到管网分叉、环网合流就开始凭感觉估算。实际工程里一个城市配气管网动辄几百个节点，手算完全不现实，这时候潮流计算就是刚需。

自己写这个6节点程序，最大的好处不是省下买软件的钱，而是能真正理解气网计算里的难点在哪。电网里电压和功率的关系已经够非线性了，天然气管道里的流量其实是压力平方差的函数，非线性程度更夸张，而且管道方向是可能反流的——一个节点周围几条管子，哪条进气哪条出气，不迭代根本不知道。

6节点是我认为最适合入门自写的规模。节点太少看不出网络效应，节点太多又会被数据准备和收敛调试拖垮。6节点可以布置成环形结构，既有主干输气，又有分支负荷，还有一种很经典的现象：环网中局部压力差导致气体从“下游”往“上游”回流。这种结果用简单手算根本发现不了，但程序一跑马上现形。

这个程序适用的场景包括：

• 天然气课程设计、毕业设计里的管网仿真部分；
• 综合能源系统研究里需要搭一个“能算气网”的基准模型；
• 工程上做小规模供气管网的方案比选，提前看压力分布是否满足用户最低供气压力；
• 作为学习牛拉法、数值雅可比、非线性方程组求解这些MATLAB基本功的载体。

说白了，这一套程序写完之后，网格加大、加入压缩机、改成动态模型，都是在现有框架上加模块，底子就是这6节点算例。

2. 算例网络与天然气潮流计算的数学表达

2.1 网络拓扑与数据准备

算例设计成一个环形配气网络，节点1是气源点，对应天然气门站出口或者调压站汇管，压力由上游调节阀控制在5.0 MPa；节点2到节点6是用气负荷节点。管道拓扑就是一个单环加弦的结构：1—2—3—4—5—6—1，其中4—5之间也直接相连，等效成了一个双环耦合的小型环网。

这种结构比纯放射状管网有意思，因为节点4周围既有来自节点3的来气，又有可能从节点5方向进气的通道，实际流向完全由压力分布决定，正好用来检验潮流计算的正确性。

节点负荷数据如下表所示。负荷单位采用工程上常用的 (10^4,\text{Nm}^3/\text{d})（每天万标准立方米），压力单位是MPa。

2.2 管道流量方程：为什么用压力平方差

天然气管道的稳态流量计算，工程上最常用的是Weymouth公式的简化形式。对于管道 (i \to j)，标况体积流量可以写成：

[ Q_{ij}=k_{ij}\cdot \mathrm{sign}(P_i-P_j)\cdot \sqrt{\left|P_i^2-P_j^2\right|} ]

其中 (k_{ij}) 是管道的导通系数，取决于管径、管长、气体相对密度、温度和压缩因子。这个方程里最关键的是“压力平方差”而不是“压力差”，因为高压气体在管道内流动时，动能项和摩擦项共同作用，最终会导出 (P^2) 的线性关系。这也是天然气潮流和电力潮流一个本质区别：电网里有功功率大致和相角差线性相关，气压网里流量和压力是平方根关系，非线性更强。

实际工程中 (k_{ij}) 通常这样算：

[ k_{ij}=C\cdot \frac{D^{2.5}}{\sqrt{\lambda L S T Z}} ]

我不建议初学者把精力耗在这个常数的量纲换算上，程序里我直接把每条管道的 (k) 值作为一个输入参数。下面的管道参数表是我按常见管网数据整理的，保证算出来的压力分布合理，量级也符合实际工程感觉。

注意，这里的 (k) 是在 (P) 单位为MPa、(Q) 单位为 (10^4,\text{Nm}^3/\text{d}) 这套单位制下的数值。如果你换成别的单位制，(k) 会差很多，这一点在第四章专门说。

2.3 节点流量平衡方程

对于任意节点 (j)，稳态下满足质量守恒：

[ \sum_{\text{流入}j} Q_i-\sum_{\text{流出}j} Q_i=L_j ]

即流入该节点的流量减去流出该节点的流量等于该节点的用气负荷。这里负荷取正值表示气体从管网被取走。

未知量怎么数？节点1压力已知，节点2到节点6共5个未知压力；节点2到节点6共5个独立的流量平衡方程。方程个数等于未知量个数，可以用牛顿-拉夫逊法求解。平衡节点1的供气量不用事先指定，等所有节点压力算出来之后，由全网流量平衡自动得到，也可以显式地作为残差补充进方程组。

3. MATLAB主程序与牛顿-拉夫逊求解框架

3.1 程序文件划分

我习惯把一个算例拆成四个文件，逻辑清晰，以后扩展也方便：

• gas6_data.m：定义节点、管道、负荷、初始值和收敛参数；
• gas6_residual.m：给定一组节点压力，计算每个负荷节点的流量平衡残差；
• gas6_jacobian.m：用数值差分计算雅可比矩阵；
• gas6_main.m：主程序，组装方程并迭代求解。

数据文件里最核心的部分是管道定义。我用的方式是每行一条管道，格式为“管道编号、起点、终点、k值”：

% gas6_data.m P1 = 5.0; % 气源/平衡节点压力 (MPa) P_init = 4.8; % 负荷节点压力初值 (MPa) loads = [0; 60; 50; 80; 70; 40]; % 节点1~6负荷 (10^4 Nm^3/d) pipes = [ 1 2 120; 2 3 95; 3 4 80; 4 5 85; 5 6 70; 6 1 105 ]; nNode = 6; nFree = nNode - 1; % 节点1为平衡节点

3.2 残差函数

残差函数是整套程序的心脏。输入是5个自由节点的压力向量，输出是一个5维残差向量，每个残差对应一个负荷节点的“流入减流出减负荷”：

function F = gas6_residual(Pfree, data) P = [data.P1; Pfree(:)]; % 组装全节点压力 loads = data.loads; pipes = data.pipes; F = zeros(length(Pfree), 1); for k = 1:size(pipes,1) a = pipes(k,1); b = pipes(k,2); kk = pipes(k,3); Q(k) = kk * sign(P(a) - P(b)) * sqrt(abs(P(a)^2 - P(b)^2)); end for j = 2:data.nNode % 节点2~6 inflow = 0; outflow = 0; for k = 1:size(pipes,1) if pipes(k,1) == j outflow = outflow + Q(k); % 以j为起点的管道，流量流出 end if pipes(k,2) == j inflow = inflow + Q(k); % 以j为终点的管道，流量流入 end end F(j-1) = inflow - outflow - loads(j); end end

这段代码没有做任何矩阵稀疏化处理，6节点规模完全够用。等以后扩展到上百节点时，再改用关联矩阵或者邻接表来组装，思路是一样的。

3.3 牛顿-拉夫逊主循环

电网潮流里牛拉法用的雅可比矩阵有明确的物理意义，气网计算可以完全照搬。方程组 (\mathbf{F}(\mathbf{P})=0) 的牛顿迭代格式为：

[ \mathbf{P}^{(k+1)}=\mathbf{P}^{(k)}-\mathbf{J}^{-1}\mathbf{F}(\mathbf{P}^{(k)}) ]

MATLAB里不建议显式求逆，直接用左除解线性方程组更快也更稳：

% gas6_main.m data = gas6_data(); Pfree = data.P_init * ones(5,1); for iter = 1:50 F = gas6_residual(Pfree, data); J = gas6_jacobian(Pfree, data); dP = J \ (-F); Pfree = Pfree + dP; if max(abs(dP)) < 1e-8 break; end end P = [data.P1; Pfree]; disp(P);

数值雅可比是我刻意选的，没有写解析导数。原因有两个：一是这个规模下数值差分速度完全不是瓶颈；二是以后改管道方程、加压缩机支路，解析雅可比要跟着推倒重来，数值雅可比只要残差函数写对，雅可比自动跟着变，省心很多。

数值雅可比的实现也很简单，用前向差分：

function J = gas6_jacobian(Pfree, data) h = 1e-6; F0 = gas6_residual(Pfree, data); J = zeros(length(Pfree), length(Pfree)); for j = 1:length(Pfree) Ppert = Pfree; Ppert(j) = Ppert(j) + h; J(:, j) = (gas6_residual(Ppert, data) - F0) / h; end end

这里有一个细节：扰动步长 (h) 不能太大也不能太小。取 (10^{-5}\sim10^{-7}) 这个区间在MPa单位制下比较稳，太小会因为浮点误差导致雅可比矩阵噪声过大，牛顿迭代反而不收敛。

4. 三个容易翻车的实现细节：方向、初值与单位

4.1 管道方向与平方根函数的“不可导尖点”

写流量计算函数时，一定要用sign(P(a)-P(b)) * sqrt(abs(P(a)^2 - P(b)^2))，而不能直接写sqrt(P(a)^2 - P(b)^2)。原因有二：第一，迭代过程中某条管道两端压力可能短暂出现 (P_a<P_b)，平方差变负，直接开根号直接NaN；第二，真实管网里管道确实可能反向输气，方向符号必须由压力差决定。

这个写法在 (P_a=P_b) 处存在不可导尖点，数值雅可比在那里会有一定误差。实际测试下来，只要不恰好卡在零流量附近反复震荡，影响很小。如果遇到某条管道流量在0附近来回跳导致不收敛，可以把步长改为带阻尼的牛顿法，比如每次迭代只走0.5倍增量，等压力场稳定下来再恢复全步长。

4.2 初值怎么给才不容易发散

牛顿法对初值敏感，气网计算尤其明显。我把所有负荷节点压力初值都设为4.8 MPa，也就是平衡节点压力的0.96倍，收敛非常顺利。但如果初值给成1 MPa，迭代很容易把压力推到负值，平方差越算越离谱，最后直接NaN。

我在调试时还试过一种更稳的办法：先用一个很小的负荷比例（比如满负荷的10%）跑一遍，把结果作为满负荷计算的初值。这就是最简单的延续法/同伦法思路，对负荷特别重的网络很有用。实际程序里如果发现满负荷直接算不收敛，可以先loads*0.1算一遍，再loads*0.5，最后算满负荷。

4.3 单位制的影响比想象中大

天然气行业单位特别混乱，MPa、bar、kPa、Pa，(10^4\text{Nm}^3/\text{d})、(\text{m}^3/\text{h})、MMSCFD都有人用。同一个物理管道，压力单位变成Pa之后，平方差就是10的12次方量级，流量公式里的 (k) 跟着变，雅可比矩阵的条件数会非常差，直接导致数值雅可比差分失真。

我踩过这个坑之后，把所有单位统一写在程序文件头部的注释里，谁接手都不会搞错。建议你也这么做：

% 单位说明： % 压力 P：MPa % 流量 Q：10^4 Nm^3/d（每天万标准立方米） % 管道长度：km（仅用于参考，不参与计算） % k 值：上述单位制下的导通系数

如果你非要用国际单位Pa和kg/s，那管道数据的 (k) 必须重新标定，不能直接从我这套数拿过去用。

5. 运行结果与环网回流现象的验证

5.1 节点压力计算结果

程序迭代收敛后的节点压力如下：

直观上看，节点4的压力是全网络最低的，因为它负荷最大（80），而且离气源最远，中间经过了好几段管道的压力损耗。节点6虽然也是“外围”节点，但它直接有一条管道连到气源节点1，所以压力反而比节点4、节点5都高，这就是网络拓扑对压力分布最直接的影响。

5.2 管道流量与“教科书级”回流现象

各管道实际流量和方向如下表。注意这里的方向是“实际流向”，不是管道定义时的参考方向：

这里面最值得玩味的是4—5管道：它实际上是反着流的，气体从节点5流向节点4，而不是从4流向5。原因不复杂：节点5从节点6方向获得了大量来气，压力稍高于节点4，于是多出来的气顺着连接管道流进了节点4，补上了节点4大负荷造成的缺口。

这个结果如果只看网络图手算，十有八九被忽略。但你拿节点平衡去校验，一切都是自洽的。也正是这个现象说明了为什么要写程序做潮流计算——环网流量分配不是按直觉走的，必须求解完整的非线性方程组。

5.3 用节点平衡校验程序正确性

判断程序算得对不对，最直接的办法是回代残差。把算出来的节点压力代回残差函数，看每个节点的“流入减流出减负荷”剩多少：

残差都在1以内，与收敛阈值 (10^{-8}) 的压力精度匹配，说明方程组确实求解到了真解附近。如果出现残差系统性地偏大，优先怀疑雅可比差分步长和单位制，其次检查管道方向处理是否写成了不合理的无符号平方根。

6. 从6节点算例向实际工程扩展的路径

6节点程序跑通之后，往上扩展其实是顺势而为。我自己常用的扩展路线是这样：

第一，节点规模扩大。把gas6_residual.m里的管道遍历改成稀疏组装，用sparse矩阵存放节点-管道关联关系，牛顿迭代里的线性求解会自动快很多。上百节点的配气管网，数值雅可比依然可以接受，前提是初值别给得太离谱。

第二，加入压缩机。压缩机是气网里的“有源元件”，类似电网里的变压器/调压器。它能把低压气体重新加压，管网方程里会出现额外的功率-流量-压比耦合关系，平衡节点也不再只有一个。本质上是在残差方程里增加一组压缩机支路方程，牛顿法框架完全不需要推倒重来。

第三，从稳态到动态。天然气管道有“管存效应”——用气低峰时管道内压力升高存储气体，高峰时压力下降释放气体，这和电网里几乎没有储能的情况完全不同。动态仿真可以把每条管道用偏微分方程离散成常微分方程组，然后用ode15s求解，程序框架上依然复用现在的节点压力向量和管道流量计算逻辑。

第四，往综合能源系统方向做。这个6节点气网模型可以和IEEE标准电网模型耦合，做一个电-气综合能源系统的联合潮流。研究里常见的“气转电”“电转气”环节，本质就是在节点负荷方程里额外增加一份耦合量。MATLAB里现在很多开源的IES工具箱，底层用的就是类似我这份程序的牛顿法内核，只是矩阵规模更大、耦合项更多。

最后说一个我自己的习惯：无论算例规模大小，我都会把节点压力、管道流量、节点残差三张表打印出来放进结果文件夹。调试时如果改了负荷或者管网拓扑，直接对比这三张表，哪里不对一眼就能看出来。这套6节点程序虽然简单，但它把天然气潮流计算的完整链路——建模、离散、求解、验证、扩展——都走了一遍，以后遇到再复杂的管网，心里也有底。

本文还有配套的精品资源，点击获取