天然气潮流计算matlab程序
1 天然气潮流计算理论
由于天然气涉及到流体的运动方程，直接计算非常复杂，因此需要提前做出一些假设来简化计算，经过研究，适当的假设对结果影响很小，因此本文对天然气系统做出如下假设：

1）管道内的天然气为等温气体，气体温度不会随着传输距离发生改变；

2）各个管道对天然气的摩擦系数相同；

（3）各个管道的气体压缩因子相同。

如果将天然气类比于电力网络潮流，则同样可以得到近似的流量方程，如式

此外，由于天然气在管道内传输时会产生压力损耗，所以每隔一段都要设置压缩机来提升节点压力，保证天然气系统的正常运行。压缩机所消耗的电能可以由式求出


然气网络中的各部分都可以用流量进行描述了。类比于电力网络，同样可以得到节点流量平衡方程，如式

管道或者节点的约束，例如流量、节点压力、压缩机消耗电能和压缩比等，如下面四式所示。

计算步骤如下：
（1）由用户负荷求得微型燃气轮机的功率𝑃g，求出其消耗的天然气流量；
（2）根据天然气节点压力与节点间流量的关系式求出前一个节点的节点压力；
（3）如果节点中存在压缩机，则可以求得压缩机消耗的等效天然气流量，将其作为一个负荷加入到潮流计算中；
（4）求出节点的节点流量方程；
（5）重复步骤（2）-（4），直至求解出整个天然气网络潮流。

2 天然气潮流计算算例结构

3 matlab程序

1）主函数

% function:计算天然气潮流
%节点类型 	标号
%普通节点 	   2,3,4,5
%恒压节点 	   1
%恒压恒流节点  6
%能计算给定各种基础数据的天然气潮流，包含各节点压力与流量

clear %清除工作空间变量
clc %清屏

%% 数据输入（标幺值）
%天然气管道参数：管道编号 长度 内径 管道常数K
Pline=[1, 500, 0.10, 0;
	   2, 400, 0.07, 0;
	   3, 300, 0.05, 0;
       4, 400, 0.07, 0;
	   5, 600, 0.07, 0];

%存储各节点流量与压力
fp=[0,   0,   0,  50,   0,   0;
    5,   0,   0,   0,   0, 5.5];

LHV=9.7;  %天然气低热值
PeMT=50;  %微型燃气轮机电功率
yitaE=0.3;%微型燃气轮机发电效率
Bk=3.8e7; %压缩机消耗等效电能算式中的常数
Zk=0.92;  %压缩机压缩因子   
   
%计算管道常数K
n=size(Pline,1);
for i=1:size(Pline,1)
    Pline(i,4)=389640*Pline(i,3)^(4/3)/sqrt(Pline(i,2)*293*0.6*0.92);
end
disp(Pline)

%计算微型燃气轮机所消耗的燃气量
fp(1,6)=PeMT/yitaE/LHV;

%% 求解天然气潮流
%主循环
for p=3:1e-5:5.5
    Hk=Bk*fp(1,6)*((fp(2,6)/p)^(Zk*0.3/1.3)-1);
    f6=1.07e-7*Hk; 
    fp(1,5)=f6+fp(1,6);
    fp(2,5)=Pi(fp(1,5),p,Pline(5,4),1);%求解出5节点压力值
    fp(2,3)=Pi(fp(1,5),fp(2,5),Pline(4,4),1);%求解出3节点压力值
    fp(2,4)=Pi(fp(1,4),fp(2,3),Pline(3,4),2);%求解出4节点压力值
    fp(1,3)=fp(1,4)+fp(1,5);%3节点总流量等于4、5节点流量之和
    fp(2,2)=Pi(fp(1,3),fp(2,3),Pline(2,4),1);%求解出2节点压力值
    fp(1,2)=fp(1,3)+15; %2节点流量为3节点流量+2节点负荷
    fp(1,1)=Pi(fp(1,2),fp(2,2),Pline(1,4),1);%求解出1节点压力值
    
    wucha=5-fp(1,1);
    if (wucha<1e-5 )
		break;
    end
end
%%  输出潮流计算结果
disp(fp);
  
。。。。。略