简介面向油气储运与能源动力专业的学生这份基于MATLAB的6节点天然气潮流计算程序可用于理解天然气输送管网中压力、流量与储存量的非线性求解过程。程序围绕状态方程、节点能量守恒、管道摩擦阻力与压降计算展开通过迭代算法如牛顿法或高斯-塞德尔法处理多方程耦合问题代码注释清晰包含主程序与管道压降函数可完整演示从建模到出结果的全链路。压缩包内共2个m文件大小仅1KB轻量精炼适合拆解研读。目前已有499人学习既可作为课堂教学的配套实验也能作为课程设计或毕业设计的前期基础帮助读者将热工流体理论灵活转化为可运行的Matlab工程代码并为进一步扩展至更大规模天然气网络分析提供参考框架。 做天然气管道水力分析的朋友大概率都遇到过这种情况上游来气压力正常管段压降也没超限但到了下游某个节点流量就是分配不均末端压力掉得离谱。这时候如果手里有一套可靠的潮流计算程序能快速算出每个节点的压力和每条管道的流量分配调试起来会省很多事。最近我把一套6节点天然气潮流计算程序用MATLAB完整实现了一遍从数学模型搭建到迭代求解再到算例验证整个过程踩了一些坑也积累了不少经验这篇就把它整理出来。这套程序解决的是天然气管网稳态运行中的核心问题给定管网拓扑、管道参数、气源供气压力和用户负荷求解各节点压力和管段流量。本质上跟电力系统里的潮流计算是对应关系。天然气系统里没有“发电厂”和“负荷”的概念但对应的是“气源点”和“用气节点”数学上都归结为非线性方程组的求解。适合正在做管网仿真、燃气输配课程设计、或者刚接触能源系统建模的工程师、研究生参考。1. 天然气管网潮流计算到底在算什么天然气潮流计算的本质是求解一组非线性方程让整个管网的流量分配和节点压力分布达到稳态平衡。跟电力系统潮流一样已知条件是管道参数、节点负荷、气源压力求解目标是管网各节点的压力和管段流量。1.1 6节点网络布局与节点分类我搭建的这套算例是一个典型的环状加枝状混合管网6个节点、7条管道。节点1是气源点连接上游天然气处理厂压力恒定控制在4.0 MPa其余5个节点均为负荷节点其中节点3、节点5、节点6接城市燃气门站节点2和节点4接工业用户。节点压力的量级需要先说一下城市高压管网入口压力一般在1.6~4.0 MPa之间工业用户要求的供气压力不低于0.8 MPa所以计算完成后要重点校核所有负荷节点的压力是否满足最低供气压力要求。这个约束在程序里我先不做硬性限制但结果分析时必须人工核查这也是工程上最常见的判断依据。1.2 初看是个简单问题仔细一想并不简单很多人第一次接触天然气潮流计算会觉得无非就是解几个压降方程把流量一步步推下去就行了。但实际做起来很快就会发现两个麻烦。第一个麻烦是环状管网的存在。如果全是枝状管网确实可以从气源点逐段向外推每个节点的压力顺着管道往下算就行。但一旦管网形成环流量分配就不再是唯一确定的同一段管道的流量可能来自两个方向必须联立求解整个方程组否则算出来就是错的。第二个麻烦是管道压降方程的高度非线性。天然气在管道里的流动压降与流量之间不是线性关系而跟流量平方、管道内径的5次方成反比关系再加上压缩因子、温度修正整个方程呈强非线性。用线性化的思路去解要么不收敛要么收敛到错误的解。所以这个计算程序的核心就是围绕一个非线性方程组用牛顿-拉夫逊迭代法在MATLAB里求解。下面逐个模块展开说。2. 数学模型把物理管网变成可求解的方程组进行潮流计算之前先要把物理管网模型化。这一步很关键很多编程上的麻烦都是模型没理清导致的。2.1 节点流量平衡方程根据质量守恒定律流入节点的天然气流量之和等于流出节点的流量之和再加上该节点的负荷。这一条和电力系统的基尔霍夫电流定律完全对应。对于每个负荷节点可以写出如下方程F_i sum(Q_in) - sum(Q_out) - L_i 0其中F_i是节点残差Q_in是从其它节点流入该节点的管道流量Q_out是流出该节点的管道流量L_i是节点负荷。对于气源节点因为压力给定它的负荷值或者说供气量是未知的由计算结果自然得到所以气源节点不列流量平衡方程只作为压力边界条件参与计算。2.2 管道压降方程与韦茅斯公式管道压降方程是天然气潮流计算的核心它的形式决定了解法的收敛性和精度。工程上常用韦茅斯Weymouth公式来近似模拟高压天然气管道的稳态流动。公式中变量比较多建议一次性定义清楚把单位统一好再写进程序。% 韦茅斯公式Q C * (P1^2 - P2^2)^0.5 % C K * (T0 / P0) * (D^2.5 / sqrt(G * T * L * Z)) % 其中 % Q - 标准状态下的体积流量m^3/d % P1 - 管道起点压力kPa % P2 - 管道终点压力kPa % D - 管道内径mm % G - 天然气相对密度空气1.0通常在0.58~0.62之间 % T - 管道内气体温度K % L - 管道长度km % Z - 压缩因子无量纲 % K - 常数取决于单位制韦茅斯公式的特點是它直接用压力平方差来表达流量避免了气体密度随压力变化带来的积分困难在中高压管网的工程计算中误差可以接受。缺点是它假设管道内流动处于完全紊流区摩擦阻力系数为常数对于低压配气管网或大管径低流速工况会有偏差。由于韦茅斯公式中流量与压力平方差是开方关系直接用它来求节点压力时会得到一个从管道方程导出的压降残差函数。在每个管道上根据两端节点压力可以计算理论流量再用理论流量与相邻管道流量差异反推节点压力或者反过来假设节点压力已知、计算每条管道流量然后检查每个节点的流量是否平衡。2.3 为什么选牛顿-拉夫逊法天然气潮流方程的求解方法有好几种常见的有方法原理优点缺点逐段推算法从气源节点逐管向外推简单直观无需迭代只适用于枝状管网节点压力迭代法固定流量迭代修正压力编程简单内存占用小收敛速度慢易振荡牛顿-拉夫逊法对残差方程组联立求雅可比矩阵迭代修正收敛快二阶收敛速度适用于环状管网需要计算偏导数编程复杂度高我最终选择牛顿-拉夫逊法理由是这套6节点管网中包含闭环结构方程组无法通过单向递推求解而牛顿-拉夫逊法在处理此类联立非线性方程组时收敛速度和稳定性都有保障。对于6节点这种小规模网络雅可比矩阵只有4到5阶计算量完全可以忽略用最直接的方式写程序反而最清晰。3. MATLAB实现从零搭一个6节点计算程序模型建立清楚后编程工作就变成了一件相对机械的事情。下面按程序模块逐一讲解我实现时的设计思路和关键代码。3.1 数据结构设计MATLAB编程里用结构体数组管理管网数据是最合理的方案。前期把每条管道的起点、终点、长度、内径、流量方向初值都定义好后面迭代求解和结果输出都可以直接索引。% 节点编号1为平衡节点气源2~6为负荷节点 % 管道数据[起点 终点 长度(km) 内径(mm) 粗糙度(mm)] pipe [ 1, 2, 12.5, 400, 0.03; 2, 3, 18.0, 350, 0.03; 3, 4, 10.0, 300, 0.03; 4, 5, 15.5, 300, 0.03; 5, 6, 8.0, 250, 0.03; 2, 5, 22.0, 300, 0.03; 3, 6, 16.0, 280, 0.03 ]; % 节点负荷[节点号 负荷(10^4 m^3/d)] load_node [ 2, 30; 3, 45; 4, 25; 5, 35; 6, 20 ];我做程序时习惯把所有的物理参数跟拓扑参数分开管理管道参数里只放几何尺寸和连接关系负荷数据单独放这样以后算不同负荷工况时不需要修改核心求解代码。3.2 节点分类与未知量提取节点分类是潮流计算里绕不开的一步。参照电力系统潮流里的节点分类法平衡节点slack node气源节点压力已知注入流量待求。对应本例的节点1。负荷节点PQ节点压力未知负荷已知。对应本例的节点2~6。对于负荷节点未知量是节点压力需要联立流量平衡方程求解。由于整个方程组里平衡节点的压力是已知的所以虽然管网有6个节点但真正需要联立求解的未知压力只有5个。程序里需要用一个索引变量把求解变量和节点号对应起来% 需计算压力的节点编号 unknown_nodes [2, 3, 4, 5, 6]; n_unknown length(unknown_nodes); % 初始压力猜测值统一给2.5 MPa P ones(6, 1) * 2500; % kPa P(1) 4000; % 平衡节点压力初值给统一值虽然不精确但对于牛顿-拉夫逊法而言只要初值不在零点附近导致雅可比矩阵奇异通常都能收敛到正确解。后面会讲初值选择对收敛性的影响。3.3 残差函数与雅可比矩阵的构建核心求解循环中每次迭代需要完成三件事根据当前节点压力计算各管道流量再计算节点流量不平衡量最后解线性方程组获得修正量。管道流量的计算函数我单独封装成了一个子函数便于复用和调试function Q pipe_flow(P_start, P_end, D, L, G, T, Z) % 韦茅斯公式计算管道流量 % 输入压力单位kPa管径mm长度km % 输出流量单位 10^4 m^3/d K 3.745e-2; % 单位换算常数不同教材取值略有差异 C K * (D^2.5) / sqrt(G * T * L * Z); Q C * sqrt(abs(P_start^2 - P_end^2)); if P_start P_end Q -Q; % 流量方向与假定方向相反 end end注意这里做了一个符号处理当终点压力大于起点压力时流量返回负值代表实际流动方向与预设方向相反。这很重要。管网的初始拓扑数据里每条管道的起点终点是我人为指定的但迭代计算过程中由于压力变化的扰动某些管道可能出现反向流如果忽略了这一点流量平衡方程会算错。在残差函数里对每个负荷节点执行流量求和function F residual(P, pipe, load_node, params) n_pipe size(pipe, 1); n_load size(load_node, 1); F zeros(n_load, 1); % 累加流入流出节点流量 for i 1:n_load node load_node(i, 1); net_flow 0; for j 1:n_pipe if pipe(j, 1) node net_flow net_flow pipe_flow(P(pipe(j,1)), P(pipe(j,2)), ... pipe(j,4), pipe(j,3), params.G, params.T, params.Z); elseif pipe(j, 2) node net_flow net_flow - pipe_flow(P(pipe(j,1)), P(pipe(j,2)), ... pipe(j,4), pipe(j,3), params.G, params.T, params.Z); end end F(i) net_flow - load_node(i, 2); end end如果说残差函数是这套程序的灵魂那雅可比矩阵就是让迭代能够高效收敛的引擎。雅可比矩阵的元素是各残差对各节点压力的偏导数。在这个程序里偏导数不手推公式而采用数值差分法计算简化编程量并保持通用性。function J jacobian_numeric(P, pipe, load_node, params) h 1e-6; % 差分步长 F0 residual(P, pipe, load_node, params); n length(F0); J zeros(n, n); unknown_nodes load_node(:, 1); for k 1:n P_temp P; P_temp(unknown_nodes(k)) P_temp(unknown_nodes(k)) h; F_plus residual(P_temp, pipe, load_node, params); J(:, k) (F_plus - F0) / h; end end数值差分的精度和步长选择有关步长太大偏导数误差大步长太小则可能被浮点截断误差淹没。经过试验1e-6在压力单位kPa、流量单位万方/天这个量级下表现稳定。3.4 牛顿-拉夫逊迭代主循环有了残差函数和雅可比矩阵迭代主循环就非常简洁了。每次迭代求解线性方程组J * delta -F然后更新压力向量。我加了一个迭代次数上限和残差范数判断条件max_iter 100; tol 1e-8; iter 0; F_norm inf; while (F_norm tol) (iter max_iter) iter iter 1; F residual(P, pipe, load_node, params); J jacobian_numeric(P, pipe, load_node, params); delta -J \ F; P(unknown_nodes) P(unknown_nodes) delta; F_norm norm(F, inf); fprintf(迭代次数: %d, 最大残差: %.6e\n, iter, F_norm); end迭代完成后节点2~6的压力就已经解出来了。再利用管道流量计算函数算出所有管道的流量和流向输出结果。实测下来这个6节点算例通常只需要5到7次迭代就能收敛到残差小于1e-8运行时间在零点几秒以内实用性完全没问题。3.5 初值选择与收敛性调整牛顿-拉夫逊法的收敛性和初值有很大关系。我在调试过程中发现初值给得太离谱比如所有节点统一给100 kPa程序偶尔会陷入不收敛。后来我把初值策略调整为在平衡节点压力基础上逐层递减给每个节点一个更贴近真实解的初始猜测。这个思路的来源是压力沿流动方向逐渐下降的物理规律% 初始压力按距离气源的管道数量递减 P zeros(6, 1); P(1) 4000; P(2) 3600; P(3) 3200; P(4) 2800; P(5) 3000; P(6) 2600;这种初值策略相当于给迭代过程一个“物理上合理”的起点显著提高了收敛稳定性和速度。对于更复杂的管网如果不知道如何给初值可以先用逐段推算法跑一遍枝状部分得到一个大致的压力分布再以此为初值做全网牛顿-拉夫逊迭代。另外一个重要的调整手段是阻尼因子。当迭代出现振荡时可以在修正量乘上一个小于1的阻尼系数alpha即P P alpha * delta。牺牲一点收敛速度换来的是稳定性的提升工程上很实用。alpha 0.8; % 阻尼因子根据收敛情况调整 P(unknown_nodes) P(unknown_nodes) alpha * delta;4. 样本结果与常见问题排查程序写成后我用一组实际参数跑了算例。这里的关键参数天然气相对密度0.6管道平均温度293 K平均压缩因子0.92。负荷取上面load_node中的数据。表里给的是我收敛后得到的主要结果。节点压力MPa负荷10^4 m^3/d1气源4.000-23.5943033.2124542.9862553.2363563.01820管段流量结果显示气源总供气量155万方/天3045253520的合计管道2-4、3-6两条管道输送量较大对应的是末端负荷较高的区域。值得注意的是节点5和节点6的压力虽然低于节点2和节点3但都高于城市门站的最低进站压力1.6 MPa要求说明这套配置下管网运行是安全的。4.1 常见问题速查表写程序的过程中我遇到过很多次计算结果明显不合理或者干脆不收敛的情况。下面这张表是踩坑经验的汇总直接列出来供参考。现象可能原因解决办法迭代发散残差越来越大初值偏离真实解太远用压力逐层递减的方法重新设置初值迭代振荡来回跳阻尼系数过大或雅可比矩阵数值差分步长不合适加入阻尼因子0.5~0.8调整差分步长某条管道流量为负压力差方向反转但公式里没加符号判断程序里加上流量方向判断负号表示反向节点压力出现负值压力平方差出现负值但开方取实数出错韦茅斯公式里用abs()后再开方并保留符号收敛到明显错误的解负荷单位不一致万方与方混用统一单位建议全部转为万方/天雅可比矩阵奇异某节点与其它节点没有管道直接相连检查拓扑数据连通性确保没有孤立节点4.2 单位制统一最容易出错的地方在天然气工程里压力单位有Pa、kPa、MPa三种常见写法流量单位有m³/d、万m³/d、Nm³/h等好几种。这套程序里我统一采用kPa作为压力单位万m³/d作为流量单位但在韦茅斯公式的常数K取值上单位换算是最容易搞错的环节。韦茅斯公式的常数跟单位制强关联不是随便抄一个就能用。标准形式的韦茅斯公式当压力单位取kPa、流量单位取m³/d时常数取0.0035左右跟气体常数、温度基准都有关系。但如果流量单位换成万m³/d常数就要相应缩万倍。我建议的做法是先把所有输入数据转换成自己设定的基准单位制常数K用与基准单位匹配的那套值并且每次换单位或换参数时就重新验算一次量纲宁可多花两分钟检查也免得后面结果全错还找不到原因。4.3 收敛判据的细节我见过不少程序直接用节点压力的修正量作为收敛判据即当两次迭代的压力差绝对值小于某个阈值时认为收敛。这种做法在这个问题上风险很大由于节点压力的数值在几千的量级如果阈值取1e-6修正量往往在初期就小于这个值但实际上流量平衡方程的残差还很大导致“伪收敛”。更好的做法是像我程序里这样用残差的无穷范数作为收敛判断标准。如果只想用压力修正量做判据阈值至少要取到1e-3以下最好同时增加流量校核条件。严格来说作为工程计算最好同时检查两种判据双条件满足才认为收敛。5. 向真实工程场景扩展的几点思考6节点只是入门。写完之后我就在想这套程序怎么做扩展才能应用到实际工程里。第一步是增加压缩机站模型。长距离输气管道每隔一段距离就需要加压站升压压缩机模型本质上是一个两端压力比与流量相关的元件加入这种元件后方程组的非线性程度会明显上升但求解框架不用变。第二步是动态仿真。稳态计算假设管网内各点参数不随时间变化但实际运行中用户负荷是24小时波动的日内调峰要求管网压力跟随负荷变化。动态仿真是在稳态模型基础上增加管道储气项、压缩机调节响应求解方法从非线性代数方程升级为微分代数方程组的时域积分。这属于进阶内容6节点模型可以作为动态仿真的空间离散初值。第三步是参数辨识与误差修正。管网运行一段时间后管道内壁粗糙度会因杂质沉积而增大实际压降会比设计值大导致原模型计算失真。利用实际运行数据反算粗糙度或压缩因子本质上是一个优化问题但底层的潮流计算程序就是目标函数的核心计算模块。这些扩展方向做完之后一套完整的天然气管网仿真平台雏形就有了。但不管扩展到哪里稳态潮流计算这个基础能力永远是最核心的积木把这套6节点程序彻底吃透后面做任何复杂模型心态都会稳很多。本文还有配套的精品资源点击获取 基于MATLAB的6节点天然气管网潮流计算程序实现与解析
简介面向油气储运与能源动力专业的学生这份基于MATLAB的6节点天然气潮流计算程序可用于理解天然气输送管网中压力、流量与储存量的非线性求解过程。程序围绕状态方程、节点能量守恒、管道摩擦阻力与压降计算展开通过迭代算法如牛顿法或高斯-塞德尔法处理多方程耦合问题代码注释清晰包含主程序与管道压降函数可完整演示从建模到出结果的全链路。压缩包内共2个m文件大小仅1KB轻量精炼适合拆解研读。目前已有499人学习既可作为课堂教学的配套实验也能作为课程设计或毕业设计的前期基础帮助读者将热工流体理论灵活转化为可运行的Matlab工程代码并为进一步扩展至更大规模天然气网络分析提供参考框架。 做天然气管道水力分析的朋友大概率都遇到过这种情况上游来气压力正常管段压降也没超限但到了下游某个节点流量就是分配不均末端压力掉得离谱。这时候如果手里有一套可靠的潮流计算程序能快速算出每个节点的压力和每条管道的流量分配调试起来会省很多事。最近我把一套6节点天然气潮流计算程序用MATLAB完整实现了一遍从数学模型搭建到迭代求解再到算例验证整个过程踩了一些坑也积累了不少经验这篇就把它整理出来。这套程序解决的是天然气管网稳态运行中的核心问题给定管网拓扑、管道参数、气源供气压力和用户负荷求解各节点压力和管段流量。本质上跟电力系统里的潮流计算是对应关系。天然气系统里没有“发电厂”和“负荷”的概念但对应的是“气源点”和“用气节点”数学上都归结为非线性方程组的求解。适合正在做管网仿真、燃气输配课程设计、或者刚接触能源系统建模的工程师、研究生参考。1. 天然气管网潮流计算到底在算什么天然气潮流计算的本质是求解一组非线性方程让整个管网的流量分配和节点压力分布达到稳态平衡。跟电力系统潮流一样已知条件是管道参数、节点负荷、气源压力求解目标是管网各节点的压力和管段流量。1.1 6节点网络布局与节点分类我搭建的这套算例是一个典型的环状加枝状混合管网6个节点、7条管道。节点1是气源点连接上游天然气处理厂压力恒定控制在4.0 MPa其余5个节点均为负荷节点其中节点3、节点5、节点6接城市燃气门站节点2和节点4接工业用户。节点压力的量级需要先说一下城市高压管网入口压力一般在1.6~4.0 MPa之间工业用户要求的供气压力不低于0.8 MPa所以计算完成后要重点校核所有负荷节点的压力是否满足最低供气压力要求。这个约束在程序里我先不做硬性限制但结果分析时必须人工核查这也是工程上最常见的判断依据。1.2 初看是个简单问题仔细一想并不简单很多人第一次接触天然气潮流计算会觉得无非就是解几个压降方程把流量一步步推下去就行了。但实际做起来很快就会发现两个麻烦。第一个麻烦是环状管网的存在。如果全是枝状管网确实可以从气源点逐段向外推每个节点的压力顺着管道往下算就行。但一旦管网形成环流量分配就不再是唯一确定的同一段管道的流量可能来自两个方向必须联立求解整个方程组否则算出来就是错的。第二个麻烦是管道压降方程的高度非线性。天然气在管道里的流动压降与流量之间不是线性关系而跟流量平方、管道内径的5次方成反比关系再加上压缩因子、温度修正整个方程呈强非线性。用线性化的思路去解要么不收敛要么收敛到错误的解。所以这个计算程序的核心就是围绕一个非线性方程组用牛顿-拉夫逊迭代法在MATLAB里求解。下面逐个模块展开说。2. 数学模型把物理管网变成可求解的方程组进行潮流计算之前先要把物理管网模型化。这一步很关键很多编程上的麻烦都是模型没理清导致的。2.1 节点流量平衡方程根据质量守恒定律流入节点的天然气流量之和等于流出节点的流量之和再加上该节点的负荷。这一条和电力系统的基尔霍夫电流定律完全对应。对于每个负荷节点可以写出如下方程F_i sum(Q_in) - sum(Q_out) - L_i 0其中F_i是节点残差Q_in是从其它节点流入该节点的管道流量Q_out是流出该节点的管道流量L_i是节点负荷。对于气源节点因为压力给定它的负荷值或者说供气量是未知的由计算结果自然得到所以气源节点不列流量平衡方程只作为压力边界条件参与计算。2.2 管道压降方程与韦茅斯公式管道压降方程是天然气潮流计算的核心它的形式决定了解法的收敛性和精度。工程上常用韦茅斯Weymouth公式来近似模拟高压天然气管道的稳态流动。公式中变量比较多建议一次性定义清楚把单位统一好再写进程序。% 韦茅斯公式Q C * (P1^2 - P2^2)^0.5 % C K * (T0 / P0) * (D^2.5 / sqrt(G * T * L * Z)) % 其中 % Q - 标准状态下的体积流量m^3/d % P1 - 管道起点压力kPa % P2 - 管道终点压力kPa % D - 管道内径mm % G - 天然气相对密度空气1.0通常在0.58~0.62之间 % T - 管道内气体温度K % L - 管道长度km % Z - 压缩因子无量纲 % K - 常数取决于单位制韦茅斯公式的特點是它直接用压力平方差来表达流量避免了气体密度随压力变化带来的积分困难在中高压管网的工程计算中误差可以接受。缺点是它假设管道内流动处于完全紊流区摩擦阻力系数为常数对于低压配气管网或大管径低流速工况会有偏差。由于韦茅斯公式中流量与压力平方差是开方关系直接用它来求节点压力时会得到一个从管道方程导出的压降残差函数。在每个管道上根据两端节点压力可以计算理论流量再用理论流量与相邻管道流量差异反推节点压力或者反过来假设节点压力已知、计算每条管道流量然后检查每个节点的流量是否平衡。2.3 为什么选牛顿-拉夫逊法天然气潮流方程的求解方法有好几种常见的有方法原理优点缺点逐段推算法从气源节点逐管向外推简单直观无需迭代只适用于枝状管网节点压力迭代法固定流量迭代修正压力编程简单内存占用小收敛速度慢易振荡牛顿-拉夫逊法对残差方程组联立求雅可比矩阵迭代修正收敛快二阶收敛速度适用于环状管网需要计算偏导数编程复杂度高我最终选择牛顿-拉夫逊法理由是这套6节点管网中包含闭环结构方程组无法通过单向递推求解而牛顿-拉夫逊法在处理此类联立非线性方程组时收敛速度和稳定性都有保障。对于6节点这种小规模网络雅可比矩阵只有4到5阶计算量完全可以忽略用最直接的方式写程序反而最清晰。3. MATLAB实现从零搭一个6节点计算程序模型建立清楚后编程工作就变成了一件相对机械的事情。下面按程序模块逐一讲解我实现时的设计思路和关键代码。3.1 数据结构设计MATLAB编程里用结构体数组管理管网数据是最合理的方案。前期把每条管道的起点、终点、长度、内径、流量方向初值都定义好后面迭代求解和结果输出都可以直接索引。% 节点编号1为平衡节点气源2~6为负荷节点 % 管道数据[起点 终点 长度(km) 内径(mm) 粗糙度(mm)] pipe [ 1, 2, 12.5, 400, 0.03; 2, 3, 18.0, 350, 0.03; 3, 4, 10.0, 300, 0.03; 4, 5, 15.5, 300, 0.03; 5, 6, 8.0, 250, 0.03; 2, 5, 22.0, 300, 0.03; 3, 6, 16.0, 280, 0.03 ]; % 节点负荷[节点号 负荷(10^4 m^3/d)] load_node [ 2, 30; 3, 45; 4, 25; 5, 35; 6, 20 ];我做程序时习惯把所有的物理参数跟拓扑参数分开管理管道参数里只放几何尺寸和连接关系负荷数据单独放这样以后算不同负荷工况时不需要修改核心求解代码。3.2 节点分类与未知量提取节点分类是潮流计算里绕不开的一步。参照电力系统潮流里的节点分类法平衡节点slack node气源节点压力已知注入流量待求。对应本例的节点1。负荷节点PQ节点压力未知负荷已知。对应本例的节点2~6。对于负荷节点未知量是节点压力需要联立流量平衡方程求解。由于整个方程组里平衡节点的压力是已知的所以虽然管网有6个节点但真正需要联立求解的未知压力只有5个。程序里需要用一个索引变量把求解变量和节点号对应起来% 需计算压力的节点编号 unknown_nodes [2, 3, 4, 5, 6]; n_unknown length(unknown_nodes); % 初始压力猜测值统一给2.5 MPa P ones(6, 1) * 2500; % kPa P(1) 4000; % 平衡节点压力初值给统一值虽然不精确但对于牛顿-拉夫逊法而言只要初值不在零点附近导致雅可比矩阵奇异通常都能收敛到正确解。后面会讲初值选择对收敛性的影响。3.3 残差函数与雅可比矩阵的构建核心求解循环中每次迭代需要完成三件事根据当前节点压力计算各管道流量再计算节点流量不平衡量最后解线性方程组获得修正量。管道流量的计算函数我单独封装成了一个子函数便于复用和调试function Q pipe_flow(P_start, P_end, D, L, G, T, Z) % 韦茅斯公式计算管道流量 % 输入压力单位kPa管径mm长度km % 输出流量单位 10^4 m^3/d K 3.745e-2; % 单位换算常数不同教材取值略有差异 C K * (D^2.5) / sqrt(G * T * L * Z); Q C * sqrt(abs(P_start^2 - P_end^2)); if P_start P_end Q -Q; % 流量方向与假定方向相反 end end注意这里做了一个符号处理当终点压力大于起点压力时流量返回负值代表实际流动方向与预设方向相反。这很重要。管网的初始拓扑数据里每条管道的起点终点是我人为指定的但迭代计算过程中由于压力变化的扰动某些管道可能出现反向流如果忽略了这一点流量平衡方程会算错。在残差函数里对每个负荷节点执行流量求和function F residual(P, pipe, load_node, params) n_pipe size(pipe, 1); n_load size(load_node, 1); F zeros(n_load, 1); % 累加流入流出节点流量 for i 1:n_load node load_node(i, 1); net_flow 0; for j 1:n_pipe if pipe(j, 1) node net_flow net_flow pipe_flow(P(pipe(j,1)), P(pipe(j,2)), ... pipe(j,4), pipe(j,3), params.G, params.T, params.Z); elseif pipe(j, 2) node net_flow net_flow - pipe_flow(P(pipe(j,1)), P(pipe(j,2)), ... pipe(j,4), pipe(j,3), params.G, params.T, params.Z); end end F(i) net_flow - load_node(i, 2); end end如果说残差函数是这套程序的灵魂那雅可比矩阵就是让迭代能够高效收敛的引擎。雅可比矩阵的元素是各残差对各节点压力的偏导数。在这个程序里偏导数不手推公式而采用数值差分法计算简化编程量并保持通用性。function J jacobian_numeric(P, pipe, load_node, params) h 1e-6; % 差分步长 F0 residual(P, pipe, load_node, params); n length(F0); J zeros(n, n); unknown_nodes load_node(:, 1); for k 1:n P_temp P; P_temp(unknown_nodes(k)) P_temp(unknown_nodes(k)) h; F_plus residual(P_temp, pipe, load_node, params); J(:, k) (F_plus - F0) / h; end end数值差分的精度和步长选择有关步长太大偏导数误差大步长太小则可能被浮点截断误差淹没。经过试验1e-6在压力单位kPa、流量单位万方/天这个量级下表现稳定。3.4 牛顿-拉夫逊迭代主循环有了残差函数和雅可比矩阵迭代主循环就非常简洁了。每次迭代求解线性方程组J * delta -F然后更新压力向量。我加了一个迭代次数上限和残差范数判断条件max_iter 100; tol 1e-8; iter 0; F_norm inf; while (F_norm tol) (iter max_iter) iter iter 1; F residual(P, pipe, load_node, params); J jacobian_numeric(P, pipe, load_node, params); delta -J \ F; P(unknown_nodes) P(unknown_nodes) delta; F_norm norm(F, inf); fprintf(迭代次数: %d, 最大残差: %.6e\n, iter, F_norm); end迭代完成后节点2~6的压力就已经解出来了。再利用管道流量计算函数算出所有管道的流量和流向输出结果。实测下来这个6节点算例通常只需要5到7次迭代就能收敛到残差小于1e-8运行时间在零点几秒以内实用性完全没问题。3.5 初值选择与收敛性调整牛顿-拉夫逊法的收敛性和初值有很大关系。我在调试过程中发现初值给得太离谱比如所有节点统一给100 kPa程序偶尔会陷入不收敛。后来我把初值策略调整为在平衡节点压力基础上逐层递减给每个节点一个更贴近真实解的初始猜测。这个思路的来源是压力沿流动方向逐渐下降的物理规律% 初始压力按距离气源的管道数量递减 P zeros(6, 1); P(1) 4000; P(2) 3600; P(3) 3200; P(4) 2800; P(5) 3000; P(6) 2600;这种初值策略相当于给迭代过程一个“物理上合理”的起点显著提高了收敛稳定性和速度。对于更复杂的管网如果不知道如何给初值可以先用逐段推算法跑一遍枝状部分得到一个大致的压力分布再以此为初值做全网牛顿-拉夫逊迭代。另外一个重要的调整手段是阻尼因子。当迭代出现振荡时可以在修正量乘上一个小于1的阻尼系数alpha即P P alpha * delta。牺牲一点收敛速度换来的是稳定性的提升工程上很实用。alpha 0.8; % 阻尼因子根据收敛情况调整 P(unknown_nodes) P(unknown_nodes) alpha * delta;4. 样本结果与常见问题排查程序写成后我用一组实际参数跑了算例。这里的关键参数天然气相对密度0.6管道平均温度293 K平均压缩因子0.92。负荷取上面load_node中的数据。表里给的是我收敛后得到的主要结果。节点压力MPa负荷10^4 m^3/d1气源4.000-23.5943033.2124542.9862553.2363563.01820管段流量结果显示气源总供气量155万方/天3045253520的合计管道2-4、3-6两条管道输送量较大对应的是末端负荷较高的区域。值得注意的是节点5和节点6的压力虽然低于节点2和节点3但都高于城市门站的最低进站压力1.6 MPa要求说明这套配置下管网运行是安全的。4.1 常见问题速查表写程序的过程中我遇到过很多次计算结果明显不合理或者干脆不收敛的情况。下面这张表是踩坑经验的汇总直接列出来供参考。现象可能原因解决办法迭代发散残差越来越大初值偏离真实解太远用压力逐层递减的方法重新设置初值迭代振荡来回跳阻尼系数过大或雅可比矩阵数值差分步长不合适加入阻尼因子0.5~0.8调整差分步长某条管道流量为负压力差方向反转但公式里没加符号判断程序里加上流量方向判断负号表示反向节点压力出现负值压力平方差出现负值但开方取实数出错韦茅斯公式里用abs()后再开方并保留符号收敛到明显错误的解负荷单位不一致万方与方混用统一单位建议全部转为万方/天雅可比矩阵奇异某节点与其它节点没有管道直接相连检查拓扑数据连通性确保没有孤立节点4.2 单位制统一最容易出错的地方在天然气工程里压力单位有Pa、kPa、MPa三种常见写法流量单位有m³/d、万m³/d、Nm³/h等好几种。这套程序里我统一采用kPa作为压力单位万m³/d作为流量单位但在韦茅斯公式的常数K取值上单位换算是最容易搞错的环节。韦茅斯公式的常数跟单位制强关联不是随便抄一个就能用。标准形式的韦茅斯公式当压力单位取kPa、流量单位取m³/d时常数取0.0035左右跟气体常数、温度基准都有关系。但如果流量单位换成万m³/d常数就要相应缩万倍。我建议的做法是先把所有输入数据转换成自己设定的基准单位制常数K用与基准单位匹配的那套值并且每次换单位或换参数时就重新验算一次量纲宁可多花两分钟检查也免得后面结果全错还找不到原因。4.3 收敛判据的细节我见过不少程序直接用节点压力的修正量作为收敛判据即当两次迭代的压力差绝对值小于某个阈值时认为收敛。这种做法在这个问题上风险很大由于节点压力的数值在几千的量级如果阈值取1e-6修正量往往在初期就小于这个值但实际上流量平衡方程的残差还很大导致“伪收敛”。更好的做法是像我程序里这样用残差的无穷范数作为收敛判断标准。如果只想用压力修正量做判据阈值至少要取到1e-3以下最好同时增加流量校核条件。严格来说作为工程计算最好同时检查两种判据双条件满足才认为收敛。5. 向真实工程场景扩展的几点思考6节点只是入门。写完之后我就在想这套程序怎么做扩展才能应用到实际工程里。第一步是增加压缩机站模型。长距离输气管道每隔一段距离就需要加压站升压压缩机模型本质上是一个两端压力比与流量相关的元件加入这种元件后方程组的非线性程度会明显上升但求解框架不用变。第二步是动态仿真。稳态计算假设管网内各点参数不随时间变化但实际运行中用户负荷是24小时波动的日内调峰要求管网压力跟随负荷变化。动态仿真是在稳态模型基础上增加管道储气项、压缩机调节响应求解方法从非线性代数方程升级为微分代数方程组的时域积分。这属于进阶内容6节点模型可以作为动态仿真的空间离散初值。第三步是参数辨识与误差修正。管网运行一段时间后管道内壁粗糙度会因杂质沉积而增大实际压降会比设计值大导致原模型计算失真。利用实际运行数据反算粗糙度或压缩因子本质上是一个优化问题但底层的潮流计算程序就是目标函数的核心计算模块。这些扩展方向做完之后一套完整的天然气管网仿真平台雏形就有了。但不管扩展到哪里稳态潮流计算这个基础能力永远是最核心的积木把这套6节点程序彻底吃透后面做任何复杂模型心态都会稳很多。本文还有配套的精品资源点击获取