简介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。节点编号类型压力初值/给定值 (MPa)负荷 ((10^4,\text{Nm}^3/\text{d}))1平衡节点/气源5.000给定由计算决定2负荷节点4.800初值603负荷节点4.800初值504负荷节点4.800初值805负荷节点4.800初值706负荷节点4.800初值402.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) 值作为一个输入参数。下面的管道参数表是我按常见管网数据整理的保证算出来的压力分布合理量级也符合实际工程感觉。管道编号起点终点参考长度 (km)参考内径 (mm)导通系数 (k)11240400120223303509533425300804452030085556202507066125350105注意这里的 (k) 是在 (P) 单位为MPa、(Q) 单位为 (10^4,\text{Nm}^3/\text{d}) 这套单位制下的数值。如果你换成别的单位制(k) 会差很多这一点在第四章专门说。2.3 节点流量平衡方程对于任意节点 (j)稳态下满足质量守恒[ \sum_{\text{流入}j} Q_i-\sum_{\text{流出}j} Q_iL_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}^{(k1)}\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_aP_b)平方差变负直接开根号直接NaN第二真实管网里管道确实可能反向输气方向符号必须由压力差决定。这个写法在 (P_aP_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) 跟着变雅可比矩阵的条件数会非常差直接导致数值雅可比差分失真。我踩过这个坑之后把所有单位统一写在程序文件头部的注释里谁接手都不会搞错。建议你也这么做% 单位说明 % 压力 PMPa % 流量 Q10^4 Nm^3/d每天万标准立方米 % 管道长度km仅用于参考不参与计算 % k 值上述单位制下的导通系数如果你非要用国际单位Pa和kg/s那管道数据的 (k) 必须重新标定不能直接从我这套数拿过去用。5. 运行结果与环网回流现象的验证5.1 节点压力计算结果程序迭代收敛后的节点压力如下节点编号压力 (MPa)15.00024.80634.67744.62654.63664.831直观上看节点4的压力是全网络最低的因为它负荷最大80而且离气源最远中间经过了好几段管道的压力损耗。节点6虽然也是“外围”节点但它直接有一条管道连到气源节点1所以压力反而比节点4、节点5都高这就是网络拓扑对压力分布最直接的影响。5.2 管道流量与“教科书级”回流现象各管道实际流量和方向如下表。注意这里的方向是“实际流向”不是管道定义时的参考方向管道实际流向流量 ((10^4,\text{Nm}^3/\text{d}))1—21 → 2165.42—32 → 3105.23—43 → 455.14—55 → 425.85—66 → 595.36—11 → 6135.2这里面最值得玩味的是4—5管道它实际上是反着流的气体从节点5流向节点4而不是从4流向5。原因不复杂节点5从节点6方向获得了大量来气压力稍高于节点4于是多出来的气顺着连接管道流进了节点4补上了节点4大负荷造成的缺口。这个结果如果只看网络图手算十有八九被忽略。但你拿节点平衡去校验一切都是自洽的。也正是这个现象说明了为什么要写程序做潮流计算——环网流量分配不是按直觉走的必须求解完整的非线性方程组。5.3 用节点平衡校验程序正确性判断程序算得对不对最直接的办法是回代残差。把算出来的节点压力代回残差函数看每个节点的“流入减流出减负荷”剩多少节点流入合计流出合计负荷残差2165.4105.2600.23105.255.1500.1455.125.80800.9595.325.87070-0.56135.295.340-0.1残差都在1以内与收敛阈值 (10^{-8}) 的压力精度匹配说明方程组确实求解到了真解附近。如果出现残差系统性地偏大优先怀疑雅可比差分步长和单位制其次检查管道方向处理是否写成了不合理的无符号平方根。6. 从6节点算例向实际工程扩展的路径6节点程序跑通之后往上扩展其实是顺势而为。我自己常用的扩展路线是这样第一节点规模扩大。把gas6_residual.m里的管道遍历改成稀疏组装用sparse矩阵存放节点-管道关联关系牛顿迭代里的线性求解会自动快很多。上百节点的配气管网数值雅可比依然可以接受前提是初值别给得太离谱。第二加入压缩机。压缩机是气网里的“有源元件”类似电网里的变压器/调压器。它能把低压气体重新加压管网方程里会出现额外的功率-流量-压比耦合关系平衡节点也不再只有一个。本质上是在残差方程里增加一组压缩机支路方程牛顿法框架完全不需要推倒重来。第三从稳态到动态。天然气管道有“管存效应”——用气低峰时管道内压力升高存储气体高峰时压力下降释放气体这和电网里几乎没有储能的情况完全不同。动态仿真可以把每条管道用偏微分方程离散成常微分方程组然后用ode15s求解程序框架上依然复用现在的节点压力向量和管道流量计算逻辑。第四往综合能源系统方向做。这个6节点气网模型可以和IEEE标准电网模型耦合做一个电-气综合能源系统的联合潮流。研究里常见的“气转电”“电转气”环节本质就是在节点负荷方程里额外增加一份耦合量。MATLAB里现在很多开源的IES工具箱底层用的就是类似我这份程序的牛顿法内核只是矩阵规模更大、耦合项更多。最后说一个我自己的习惯无论算例规模大小我都会把节点压力、管道流量、节点残差三张表打印出来放进结果文件夹。调试时如果改了负荷或者管网拓扑直接对比这三张表哪里不对一眼就能看出来。这套6节点程序虽然简单但它把天然气潮流计算的完整链路——建模、离散、求解、验证、扩展——都走了一遍以后遇到再复杂的管网心里也有底。本文还有配套的精品资源点击获取