简介一份基于MATLAB的6节点天然气潮流计算教学程序面向能源动力及相关专业初学者用于理解天然气网络压力、流量与储存量的求解方法。程序以6节点简化模型为对象完整涉及状态方程、能量守恒、管道阻力与压降计算等核心环节压降估算可参考达西-韦斯巴赫公式或海曾-威廉公式并采用牛顿法或高斯-塞德尔迭代求解相互关联的非线性方程组帮助读者把理论公式落地为可运行代码。压缩包内含2个M脚本文件体积仅约1KB代码注释清晰、结构紧凑包含网络数据输入、潮流计算主体、结果输出与可视化流程。已有499人学习下载适合作为天然气仿真分析入门的参考模板也可在此基础上进一步扩展至更复杂的管网模型。 做天然气潮流计算的朋友应该都有这种体会书上的公式看着不难真正能把一个算例跑通、跑对需要跨过的坑比想象中多得多。最近帮师弟调试课程设计正好用的是这套6节点天然气潮流计算MATLAB程序趁着记忆还热乎把这个算例从建模思路、数据准备到程序实现和调试经验完整梳理一遍。内容不绕弯子直接把能复现的算例数据、代码框架和踩坑记录放在这里给同样在做天然气管网仿真、综合能源系统潮流计算或者刚接触MATLAB编程的同学参考。1. 天然气潮流计算到底在算什么1.1 先把问题跟电力潮流放一起看很多人第一次看到“潮流计算”四个字会下意识把它理解成电力系统里的牛拉法、PQ分解法那一套。这个直觉没错天然气管网潮流和电力系统潮流本质上是一类问题给定网络拓扑、边界条件和负荷求解全网的状态变量分布。电力潮流求的是节点电压幅值和相角天然气潮流求的是节点压力和管段流量。两者的不同点也很明显。电力系统里有功潮流跟相角差近似线性无功跟电压差强相关所以可以分PQ分解天然气管网里管段流量跟压降之间的关系是非线性的平方关系而且没有“相角”这类缓变量方程组的非线性程度更高对初值更敏感。换句话说电力潮流不收敛的时候通常调调初值就好天然气潮流不收敛的时候你可能要回头检查管段方程是不是写错了。但有趣的是求解框架完全可以复用列节点守恒方程加元件特性方程组成非线性方程组后用牛顿-拉夫逊法迭代求解。所以如果你已经写过电力潮流程序天然气潮流对你来说就是个“换了元件方程”的版本。1.2 管网模型里两个核心方程天然气管网模型的本质由两类方程构成。第一类是节点流量平衡方程。对任意一个节点流入该节点的流量之和减去流出该节点的流量之和必须等于该节点的负荷用气量如果有注入源比如气源节点则注入量减去全部流出量等于零。数学上写成[ \sum_{j \in N(i)} q_{ij} Q_{load,i} ]这个方程的意义跟电力潮流的KCL一模一样——物质守恒。在程序里它是构造残差向量的基础。第二类是管段压降方程。天然气在管道里流动时由于摩擦阻力产生压降。对一段连接节点 (i) 和 (j) 的管道简化的压降-流量关系可以写成[ p_i^2 - p_j^2 R_{ij} \cdot q_{ij}^2 ]这里 (R_{ij}) 是管段的等效阻力系数跟管长、管内径、天然气物性相对密度、压缩因子、温度有关实际工程里常用Weymouth公式或Darcy公式计算。这个方程最大的特点是压力出现在平方项里而流量也出现在平方项里典型的双向非线性。正因为这个“平方对平方”的关系很多教材把状态变量直接取成压力的平方 (P2 p^2)方程形式会简洁不少迭代时也更好处理。提示潮流计算中压力一定要用绝压不能用表压。新手经常在这一点上翻车算出来的结果偏差大到离谱。1.3 为什么选6节点作为入门算例我先解释一下为什么拿6节点而不是3节点或者IEEE标准节点来说事。3节点算例只能做单一链式结构气源-中游-末端体现不出环网的流量分配问题而IEEE那种动辄几十节点的算例对初学者来说数据准备和结果校核都是负担。6节点是一个非常适中的规模结构上可以同时包含串联干线、分支支路和环网能完整反映天然气管网的主要特征结果规模又足够小可以手工校核任意一条管段的压降计算是否合理。很多教材和论文里的演示算例都采用这个规模参考资料也容易找。我自己带学生的感受是6节点算例跑通之后理解“平衡节点”“定流量节点”“迭代初值”“残差收敛”这些概念就都具体化了再往上扩到十几节点或几十节点只是数据量的增加不再是思路层面的障碍。2. 6节点算例的搭建与数据准备2.1 管网拓扑结构设计我设计的这个6节点算例拓扑上力求覆盖三种典型结构链式干线、分支支路、以及一个闭合环网。具体结构是这样的节点1是气源节点相当于电力系统中的平衡节点压力给定节点1到节点2到节点3是主输气干线管径较大节点2分出支路到节点4节点4再向下到节点5这是末端支路节点3和节点4之间通过节点6形成一个闭合环路即3-6-4-3。为什么要布置一个环网因为环网是天然气潮流计算中最容易出问题的地方流量分配不是直观的“上游到下游”而是需要解方程组才能确定流量方向也可能跟初始猜测相反。如果只用树状管网用序贯法手算也能算体现不出牛顿-拉夫逊法的价值。环网结构一旦出现就必须用到雅可比矩阵迭代求解。2.2 节点负荷与管段参数表为了让这个算例能直接复现我把数据全部列出来。单位系统统一使用压力单位MPa流量单位百万立方米每天MMscmd管道阻力系数 (R) 的单位为 (\text{MPa}^2 / (\text{MMscmd})^2)。节点负荷数据节点编号节点类型负荷MMscmd说明1平衡节点气源0压力固定为5.0 MPa2负荷节点0.8民用/工业用气3负荷节点1.2民用/工业用气4负荷节点0.5支路用户5负荷节点0.6末端用户6负荷节点0.7环网用户总负荷是3.8 MMscmd全部由节点1的气源供应。管段参数和等效阻力系数管道编号起点终点管长km管径mm等效阻力系数R11282000.18223102000.2232461501.1044581501.4553651500.9266461501.10我给的R值是折算后的等效值实际工程中需要用Weymouth公式从管长、管径、天然气物性计算。这里直接给折算值的目的是先把潮流计算的框架打通R的精确计算放到后面再完善。如果读者打算换成实际公式只需要写一个从管段参数计算R的函数然后替换掉输入数据中的R列即可。2.3 参数单位换算这个坑预算是非常关键的一步也是最容易出大问题的地方。我第一次帮别人调这个程序的时候发现计算结果跟手算完全对不上查了一个多小时最后发现是压力单位混用了迭代内部用的Pa但管道阻力系数是用MPa代入公式算的两者差了12个量级程序不跑飞才怪。建议从一开始就确定统一的单位体系并且在整个程序内保持一致。我习惯的做法是外部数据输入用工程单位MPa、km、mm进入程序后在读取数据时立刻转换成迭代内部单位Pa、m算完再转回工程单位输出。这样虽然多两步转换代码但能避免最常见的单位灾难。3. MATLAB程序核心实现3.1 输入数据的代码组织在MATLAB里我习惯用向量存节点数据、用矩阵存管段数据。不需要复杂结构体数据量小的时候直接数组操作更快、更直观。对应上面的算例输入模块长这样% 节点负荷单位MMscmd nodeLoad [0; 0.8; 1.2; 0.5; 0.6; 0.7]; % 管段数据每行 [起点, 终点, 等效阻力系数R] pipe [ 1 2 0.18 2 3 0.22 2 4 1.10 4 5 1.45 3 6 0.92 6 4 1.10 ]; % 气源节点和压力 sourceNode 1; sourcePressure 5.0; % MPa % 迭代参数 tol 1e-8; % 收敛容差 maxIter 100; % 最大迭代次数这样组织数据的优点是管段数量和拓扑修改都只动这两个数组程序其他部分完全不用改。等到以后扩展到实际管网把Excel或数据库的数据读进来替换这里的赋值语句就行。3.2 牛顿-拉夫逊法迭代框架核心求解思路是把状态变量设为各节点压力平方 (P2 p^2)。对每个节点气源节点除外构造残差[ F_i \sum_{j \in N(i)} q_{ij} - Q_{load,i} ]其中管段流量由压降方程反解[ q_{ij} \text{sign}(P2_i - P2_j) \cdot \sqrt{\frac{|P2_i - P2_j|}{R_{ij}}} ]加 (\text{sign}) 的原因是环网中流量方向未知必须根据两端压力大小自动判断方向。雅可比矩阵的元素是 (F_i) 对 (P2_i) 的偏导数。对一段管道的流量项能直接推到解析表达式[ \frac{\partial q_{ij}}{\partial (P2_i)} \frac{1}{2 \sqrt{|P2_i - P2_j| \cdot R_{ij}}} ]在程序里组装雅可比矩阵时对每根管道这部分的贡献同时作用于四个位置(J(i,i))、(J(j,j)) 加正号(J(i,j))、(J(j,i)) 加负号。核心迭代代码的骨架如下P2 sourcePressure^2 * ones(6,1); % 初值 for iter 1:maxIter F zeros(6,1); J zeros(6,6); for k 1:size(pipe,1) i pipe(k,1); j pipe(k,2); R pipe(k,3); dp2 P2(i) - P2(j); q sign(dp2) * sqrt(abs(dp2) / R); F(i) F(i) q; F(j) F(j) - q; dq 1 / (2 * sqrt(abs(dp2) * R 1e-12)); J(i,i) J(i,i) dq; J(j,j) J(j,j) dq; J(i,j) J(i,j) - dq; J(j,i) J(j,i) - dq; end F F - nodeLoad; % 节点流量平衡残差 % 处理平衡节点将气源节点方程替换为定压约束 F(sourceNode) P2(sourceNode) - sourcePressure^2; J(sourceNode, :) 0; J(sourceNode, sourceNode) 1; % 迭代修正 delta -J \ F; P2 P2 delta; if max(abs(delta)) tol break; end end % 换算回压力 pressure sqrt(P2);这里有个细节值得注意在求导表达式中我加了 (1e-12) 作为数值保护。当某管段两端压力平方差接近零时导数会趋于无穷大加一个极小量可以防止数值溢出。虽然6节点算例里不一定触发这个问题但这个习惯在扩展到大型管网时能救命。3.3 平衡节点的处理逻辑平衡节点气源的处理是牛顿-拉夫逊法在这个问题里的关键细节。节点1的压力是给定的所以它的节点方程不再是流量平衡方程而是定压约束方程[ P2_1 p_{source}^2 ]在矩阵里做起来很直接把该节点对应的残差F行改成 (P2_1 - p_{source}^2)把雅可比矩阵对应行清零只在对角元置1。这样既保留了矩阵结构的完整性又实现了压力约束。有个常见错误是直接删掉平衡节点所在行和列这在网络规模大、稀疏矩阵处理时反而引入麻烦而且在输出节点压力时还得再映射回去代码绕来绕去容易出bug。固定行替换的方法在6节点规模下迭代次数几乎没差别代码还更简洁。3.4 收敛判据怎么选我用的收敛判据是压力平方修正量的最大绝对值小于 (10^{-8})。这个阈值在压力以MPa为单位时对应压力修正大约在 (10^{-9}) MPa量级精度完全够用。也有一些实现用残差范数做判据比如 (||F||_\infty 10^{-6})。两种判据在已经收敛的迭代序列上差别不大但用状态修正量做判据有个好处不会因为某个节点负荷特别小导致残差天然很小而产生“假收敛”。实际调试时建议两种判据都打印出来看一个判断解的稳定性一个判断方程满足程度。4. 运行结果怎么看4.1 参考输出示例按上面的数据和代码跑完压力结果大致如下节点编号压力MPa15.00024.8334.6744.7154.5564.62整体规律符合物理直觉离气源越远、负荷越重的节点压力越低环网中节点3、4、6的压力比较接近这正是环网平衡流量后自然形成的压力分布。末端节点5压力最低因为它在支路末端且管径偏小。如果算出来的结果里出现某个节点压力比气源还高那一定有问题——可能是流量方向写反了也可能是雅可比矩阵符号错误。4.2 自己动手校核的三种方法程序跑通不等于跑对我建议每个结果都要做三个层面的校核。第一是全局校核。所有负荷相加是3.8 MMscmd气源节点的供气量必须也是3.8通过节点1的流量平衡方程反算。如果不相等说明某个节点方程写错了。第二是局部校核。随便挑一根管道比如管段1-2用输出结果计算 (p_1^2 - p_2^2)再计算 (R \cdot q^2)两边应该严格相等。这是最直接的方程校验能快速定位程序里哪一步出了问题。第三是灵敏度校核。把某个节点的负荷增大10%看这个节点以及上游节点的压力是否下降。如果压力反而上升了那程序肯定有逻辑错误。这个方法不需要额外工具就是一种很好的程序自检习惯。4.3 从结果反推管网瓶颈潮流计算不只是“算个压力”而已它的输出可以直接用于管网规划与运行分析。比如这个算例的结果里节点5的压力是4.55 MPa如果设计要求末端压力不低于4.6 MPa那这个管网就是不满足要求的需要采取升压或扩容措施。这种“结果驱动决策”的思路是潮流计算作为工具的真正价值在管网规划阶段预测不同负荷场景下的压力分布在运行阶段判断管网是否有瓶颈、是否需要增压在综合能源系统研究中为电-气耦合分析提供天然气侧的运行状态。6节点算例虽然小但这个分析逻辑跟实际工程完全一致。5. 调试实录那些年踩过的坑5.1 不收敛先查初值和符号牛顿法不收敛的原因按出现频率排序大概是初值偏离太远、雅可比矩阵符号写反、单位混用。初值问题是最常见的。这个算例里气源压力5.0 MPa如果你把所有节点初始压力设为0迭代第一步就会遇到很大的导数矩阵数值特性差很容易震荡发散。我建议把初始压力设为气源压力的0.8倍左右也就是初步认为不同节点压力差别不会太大。对绝大多数中压、低压管网这个初值策略都很稳。符号问题则隐蔽得多。雅可比矩阵的符号错误会让残差序列呈现“来回蹦”的特征一会儿正、一会儿负中间还出现过接近收敛的情况然后突然又弹出很远。遇到这种症状最有效的排查办法是把迭代前几步的雅可比矩阵打印出来手工检查第一个节点的对角线元素的导数表达式是否正确。5.2 压力变成负数或复数如果你在迭代过程中发现压力平方出现负值说明某个节点的“压力平方”被迭代修成了负数开方后就是复数MATLAB会输出NaN或者复数结果。这种情况绝大多数是物理模型本身出了问题负荷过大、管径过小、气源压力不足。在6节点算例里如果你把总负荷加大到10 MMscmd以上很可能触发这个问题。它的含义是在这个负荷水平下现有管网结构根本无法满足供气要求潮流计算在数学上无解。遇到这种情况不要试图靠调初值绕过而应该回头检查参数合理性。这也是潮流程序的一个隐形价值它可以作为管网可行性的判定工具。5.3 常见问题速查表现象可能原因处理办法完全不收敛残差猛增初值太差或雅可比矩阵符号错误调初值为气源压力的0.8倍打印雅可比矩阵逐项核对迭代在某个值附近振荡流量方向误判sign函数缺失检查管段流量计算是否加了sign节点压力为NaN或复数负荷过大或管径过小无可行解减小负荷或增大管径检查R值结果跟手算对不上单位混用统一用MPa和MMscmd或迭代内部统一用SI单位气源节点流量不为总负荷平衡节点方程处理错误检查该节点的定压方程是否被错误替换掉某根管道流量为负但压力差为正管段端点编号顺序与实际方向不同流量方向以sign为准负号是正常现象这个小表基本覆盖了我自己调试这个程序时遇到的主要问题也几乎覆盖了学生在这类作业里问过我的所有问题。6. 从6节点往外扩展的方向6.1 换数据就能跑的拓展思路6节点程序跑通之后往实际管网扩展的路径其实很清晰。把节点数、管道数从数组里替换成实际数据加上压缩机模型、阀门模型再把稀疏矩阵技术用起来就是一个简化版的天然气管网仿真工具。我实际测试过把同样的代码框架扩展到20节点左右的环形管网计算时间仍然在毫秒级。真正需要改进的是矩阵求解部分当节点数到几百甚至上千时直接把雅可比矩阵存成稠密矩阵就不合适了需要改成稀疏存储并配合稀疏线性方程求解。不过这些都属于“工程优化”底层方程和求解逻辑跟6节点算例没有本质区别。所以把这套基础程序吃透等于给后续拓展打好了地基。6.2 综合能源系统方向如果做的是电-气耦合系统研究可以把这份6节点天然气潮流程序跟一个简单的电力系统潮流程序接力起来用天然气网输出的气源供气量去约束燃气轮机的出力再用电网的负荷需求去影响天然气的用气量。两步交替迭代就是一个最简单的电-气联合潮流雏形。这个方向这几年在综合能源系统研究里很热门但很多论文里的算例就是把这两个程序的接口做了一下对接。基础能力还是各自领域的潮流求解所以先练好6节点天然气潮流对后面做系统级研究会很有帮助。6.3 一点实际操作的体会最后分享一个我每次带新手都会强调的建议跑通这个程序之后别急着交差把负荷数据改一改再跑几遍。比如把节点5的负荷从0.6改成1.5观察节点压力怎么变把节点4和节点6之间的连接断开看环网变成支路后流量怎么重新分配。这些“折腾”对理解管网运行特性的帮助比单纯把程序跑通要大得多。我最初带师弟做这个题目时他调试到凌晨三点最后发现是压力单位混用。折腾的过程虽然痛苦但那次之后他对单位换算、初值选择这类细节有了刻骨铭心的记忆后面再调电网、燃气系统的程序都顺了很多。这种基本功靠看教程是学不来的必须亲自踩一次坑。本文还有配套的精品资源点击获取