一种用于综合能源网中热力管网运行参数的动态计算方法

摘要

本发明涉及一种用于综合能源网中热力管网运行参数的动态计算方法，属于综合能源网数字仿真技术领域。本发明的热力管网运行参数的动态计算方法，充分考虑了热力管网中工质物性变化、热力管网动态过程中水力工况和热力工况间的相互影响，不仅适用于热水等热力工况对水力工况影响较小的动态过程，也能适用于导热油、蒸汽等热力工况对水力工况影响较大的动态过程。同时，本发明方法从数学描述和表达形式上完成了统一，可实现解算方法的统一，便于计算机编程实现，利用矩阵进行线性求解，具有很好的应用前景。

权利要求

1.一种用于综合能源网中热力管网运行参数的动态计算方法，其特征在于该方法包括以下步骤：

设定：综合能源网中热力管网中的节点总数为m，支路总数为n，热力计算步长为Δt，水力计算步长为ΔT，水力求解系数矩阵为A、流量补偿向量为B，热力求解系数矩阵为D、热力补偿向量为E，初始化系数矩阵A、D的元素及补偿向量B、E元素为0，初始化各节点压力Pi，为

0.1、节点温度Ti为20，i＝1,2…m，根据设计工况初始化各支路流量wj、各支路与外界的热交换功率qj及各节点的海拔高度hi，j＝1,2…n；根据综合能源网中热力管网的拓扑结构,得到二维数组AI[i][ai]和BI[i][ib],其中二维数组AI[i][ai]的行号i为节点编号，第i行的行元素为与流入节点i的支路相连的相邻节点编号，ai为与流入编号为i的节点的支路相邻的节点总数，二维数组BI[i][ib]的行号i为节点编号，第i行的行元素为与流出编号为i的节点的支路相连的相邻节点编号，ib为与流出编号为i的节点的支路相邻的节点总数；

(1)根据热力管网中任意节点i内输送介质的温度Ti和介质压力Pi，i＝1,2…m，m为综合能源网中热力管网中的节点总数，利用介质不同温度对粘度的映射函数f1，根据温度、压力求介质密度的函数f2和求介质比热的函数f3，分别计算支路j中输送介质的动力粘度νj、密度ρbj以及第i号节点中的输送介质的密度ρi和比热ci；

νj＝f1(Tcj)

ρbj＝f2(Pcj,Tcj)

ρi＝f2(Pi,Ti)

ci＝f3(Pi,Ti)

cj＝f3(Pcj,Tcj)

式中，j＝1,2…n，n为综合能源网中热力管网中的支路总数，cj为与支路j相连的上节点编号；

(2)根据热力管网中任意支路j内输送介质的温度Tj，j＝1,2…n，n为综合能源网中热力管网中的支路总数，结合步骤(1)求得的热力管网中任意支路j的流体动力黏度νj，利用下式求出支路j的最大通流系数Kvj及支路j的内容积VBj；

式中，lj为支路j的管线长度，uj为支路j的管线中输送介质的流速，dj为支路j的管线的当量直径，Rej为支路j中介质流动时的流态判断变量雷诺数,λj为支路j的管线的沿程阻力系数,ρbj为支路j中介质的密度；

(3)从调门曲线选取三个调门中间开度θ1,θ2,θ3，根据调门曲线选取的与开度θ1,θ2,θ3相对应的调门通流比例χ1,χ2,χ3，计算支路j的通流比例χ和实际通流系数Cvj：对支路j进行判断，若支路j不带有调节阀门，则将支路j的通流比例χ置为1，若支路j带有调节阀门，则采用四次多项式拟合的方法，其中，aa0,aa1,aa2,aa3,aa4分别为拟合多项式系数，对阀门通流比例χ进行多项式拟合计算：

2 3 4

χ＝aa0+aa1θ+aa2θ+aa3θ+aa4θCvj＝χKvj

上式中，Kvj为步骤(2)中支路j的最大通流系数，χ为调门的铜留比例，取值为0‑1，θ为调节阀门的实际开度，θ的取值范围为0‑1；

(4)根据与支路j相对应的管线内的流体质量流量wj，利用下式，计算支路j的流导Rj：上式中，ρbj为步骤(1)的支路j中流动介质的密度；

(5)根据支路j中流体质量流量wj、与支路j相连的上节点的海拔高度hcj、与支路j相连的下节点的海拔高度hdj，利用下式，计算与热力管网水力迭代计算过程中支路j相对应的流量补偿项Cj：

式中，g为重力加速度，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；

(6)根据步骤(4)的支路j流导Rj，利用下式，计算得到一个用于热力管网水力计算的水力求解系数矩阵A的对角元素aii以及与节点i相对应的非对角元素aik,aki：aik＝aki＝Rjio

式中，二维数组AI[i][ai]的行号i为节点编号，第i行的行元素为与流入编号为i的节点的支路相连的相邻节点编号，ai为与流入编号为i的节点的支路相邻的节点总数，二维数组BI[i][ib]的行号i为节点编号，第i行的行元素为与流出编号为i的节点的支路相连的相邻节点编号，ib为与流出编号为i的节点的支路相邻的节点总数，下标CI[i][ai+ib]为一个二维数组，二维数组中的行号i为节点编号，第i行的行元素依次为ai个流入编号为i的节点的支路编号和ib个流出编号为i的节点的支路编号，iki为遍历流入编号为i的支路的计数变量，iki＝1,2…ai，iko为遍历流出支路i的计数变量，iko＝1,2…ib，jio(jio＝CI[i][1],CI[i][2]…CI[i][ai+ib])为与节点i的相连的支路编号，k为与节点i相连的编号为jio的支路相连的相邻节点编号，k＝AI[i][1],…AI[i][ai],BI[i][1],…BI[i][ib]；

(7)根据步骤(5)求得的热力管网各支路的流量补偿项Cj，利用下式，计算得到一个流量补偿向量B，流量补偿向量B中的元素为bi，i＝1,2…m：式中，iki为遍历流入编号为i的节点的支路的计数变量，iki＝1,2…ai，iko为遍历流出节点i的支路的计数变量，iko＝1,2…ib；

(8)根据步骤(6)和步骤(7)计算的水力求解系数矩阵A和流量补偿向量为B的元素值，得到一个热力管网水力计算的线性方程组，利用高斯‑赛德尔迭代的方法求解该线性方程t+ΔT

组，得到t至t+ΔT计算时段内热力管网中各节点压力Pi ，i＝1,2…m：t+ΔT

式中，Pi 为t至t+ΔT计算时段内热力管网中节点i的压力；

(9)记录上一计算时步中热力管网内各节点压力Pi，i＝1,2…m，利用下式计算在t至t+ΔT计算时段内各节点压力与t‑ΔT至t计算时段内各节点压力差ΔPi，选取所有节点压力t+ΔT

误差值中的最大值作为收敛判据ε，并将更新的各节点压力Pi 赋值给Pi：t+ΔT

ΔPi＝Pi ‑Pi

ε＝max{ΔP1,ΔP2…ΔPm}t+ΔT

Pi＝Pi

(10)根据步骤(9)的热力管网各节点的压力Pi、步骤(4)的热力管网各支路的流导Rj及步骤(5)的各支路流量补偿项Cj(j＝1,2…n)，利用下式更新支路j中的介质质量流量wj；

wj＝Rj(Pcj‑Pdj)‑Cj式中，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；

(11)设置一个热力管网的水力迭代计算阈值δ，对步骤(9)的进行判断：若ε≥δ则返回步骤(3)，若|ε|＜δ，则进入步骤(12)；

(12)根据步骤(10)的介质质量流量wj及步骤(1)中求得的各支路中流动介质的比热cj,得到一个用于热力管网的热力计算的热力求解系数矩阵D，热力求解系数矩阵D中与节点i(i＝1,2…m)相对应的非对角元素dif为：dif＝wjocjoΔt

式中，jo为流出编号为i的节点的支路编号,jo＝CI[i][ai+1]),CI[i][ai+2]…CI[i][ai+ib]，f为与流出节点i的支路jo相连的节点编号，f＝BI[i][1],…BI[i][ib]；

(13)根据步骤(2)的各支路对应管线的内容积VBj、步骤(1)的各节点处流动流体的密度ρi和比热ci，i＝1,2…m，利用步骤(10)的热力管网各支路流动介质质量流量wj，利用下式得到热力管网中各节点的当量容积Vi以及步骤(12)的用于热力管网的热力求解系数矩阵D的对角线元素dii：

式中，Vi为编号为i的节点的当量内容积，ikio为遍历与编号为i的节点相连的支路的计数变量，ikio＝1,2…ai+ib，iki为遍历流入编号为i的节点的支路的计数变量，iki＝1,2…ai：

(14)根据热力管网中各支路对外界的热功转换功率qj，j＝1,2…n，以及步骤(13)的热力管网中各节点的当量容积Vi、步骤(1)的热力管网各节点处流动介质的密度ρi和比热ci，利用下式得到热力管网各节点处注入的热功率qni：并得到一个热力补偿向量为E，热力补偿向量为E中的元素ei：ei＝‑Viρici+Δtqnj式中，iki为遍历流入编号为i的节点的支路的计数变量；

(15)根据步骤(12)‑(14)的热力计算系数矩阵D和热力补偿向量E的元素，得到用于热力管网热力计算的线性方程组，利用高斯‑赛德尔迭代方式求解该线性方程组，得到t至t+t+Δt

Δt计算时段内热力网络各节点处输送介质的温度Ti ：t+Δt

(16)根据步骤(15)的各节点流动介质的温度Ti ，i＝1,2…m，利用下式计算出热力管网热力计算的t至t+Δt计算时段内各节点中的介质温度与t‑Δt至t计算时段内各节点介质温度差ΔTi，选取所有节点介质温度误差值中的最大值作为收敛判据ε'，并将更新的各t+Δt

节点温度Ti 赋值给Ti：

t+Δt

ΔTi＝Ti ‑Ti

t+Δt

Ti＝Ti

ε'＝max{ΔT1,ΔT2…ΔTm}；

(17)设置一个热力管网的节点温度的迭代计算阈值δ'，对步骤(16)ε'，进行判断：若|ε'|≥δ'，则返回步骤(2)，若|ε'|＜δ'，则结束计算，并将步骤(9)得到的热力管网内各节点处介质的压力Pi、步骤(10)得到的热力管网内各支路中的介质质量流量wj、步骤(16)得到的各节点出流动介质的温度Ti作为综合能源网中热力管网的运行参数，实现能源网中热力管网运行参数的动态计算。

说明书

一种用于综合能源网中热力管网运行参数的动态计算方法

技术领域

[0001] 本发明涉及一种用于综合能源网中热力管网运行参数的动态计算方法，属于综合能源网数字仿真技术领域。

背景技术

[0002] 综合能源网覆盖了热力侧子网、机电侧子网、燃料侧子网等不同类型能量网络及其具有不同响应时间的设备部件等，热力管网的动态仿真技术在综合能源系统仿真中起到
越来越重要的作用。

[0003] 目前的热力管网动态仿真建模主要从流体动力学的角度完成水力工况的建模，针对复杂的热力工况、工质物性变化、相变等因素缺乏统一的处理手段，且当热力管网和电力
网联立解算的时候因建模方法和数学表达不统一无法统一其算法，程序实现困难，为综合
能源系统仿真建模带来困扰。如导热油输送管网在导热油温度上升时，未考虑到导热油的
物性影响，随着温度上升，导热油因为动力粘度的变化引起流阻的变化，热水输送过程中，
未考虑因热水温度上升可能引起相变的影响，忽略了热力工况对水力工况的影响，孤立的
对待热力工况和水力工况。

发明内容

[0004] 本发明的目的提出一种用于综合能源网中热力管网运行参数的动态计算方法，借鉴机电侧子网机‑网建模的方法，根据传统热力学理论从“广义电路”的角度对热力管网模
型进行统一描述和数学表达，完成综合能源系统动态仿真模型求解的算法统一和程序实
现。

[0005] 本发明提出的用于综合能源网中热力管网运行参数的动态计算方法，包括以下步骤：

[0006] 设定：综合能源网中热力管网中的节点总数为m，支路总数为n，热力计算步长为Δt，水力计算步长为ΔT，水力求解系数矩阵为A、流量补偿向量为B，热力求解系数矩阵为D、
热力补偿向量为E，初始化系数矩阵A、D的元素及补偿向量B、E元素为0，初始化各节点压力
Pi(i＝1,2…m)为0.1、节点温度Ti为20，根据设计工况初始化各支路流量wj(j＝1,2…n)、各
支路与外界的热交换功率qj及各节点的海拔高度hi；根据综合能源网中热力管网的拓扑结
构,得到二维数组AI[i][ai]和BI[i][ai],其中二维数组AI[i][ai]的行号i为节点编号，第
i行的行元素为与流入节点i的支路相连的相邻节点编号，ai为与流入编号为i的节点的支
路相邻的节点总数，二维数组BI[i][ai]的行号i为节点编号，第i行的行元素为与流出编号
为i的节点的支路相连的相邻节点编号，ib为与流出编号为i的节点的支路相邻的节点总
数；

[0007] (1)根据热力管网中任意节点i内输送介质的温度Ti和介质压力Pi，i＝1,2…m，m为综合能源网中热力管网中的节点总数，利用介质不同温度对粘度的映射函数f1，根据温度、
压力求介质密度的函数f2和求介质比热的函数f3，分别计算支路j中输送介质的动力粘度
νj、密度ρbj以及第i号节点中的输送介质的密度ρi和比热ci；

[0008] νj＝f1(Tcj)

[0009] ρbj＝f2(Pcj,Tcj)

[0010] ρi＝f2(Pi,Ti)

[0011] ci＝f3(Pi,Ti)

[0012] cj＝f3(Pcj,Tcj)

[0013] 式中，j＝1,2…n，n为综合能源网中热力管网中的支路总数，cj为与支路j相连的上节点编号；

[0014] (2)根据热力管网中任意支路j内输送介质的温度Tj，j＝1,2…n，n为综合能源网中热力管网中的支路总数，结合步骤(1)得到热力管网中任意支路j的流体动力黏度νj，利
用下式求出支路j的最大通流系数Kvj及支路j的内容积VBj；

[0015]

[0016]

[0017]

[0018]

[0019] 式中，lj为支路j的管线长度，uj为支路j的管线中输送介质的流速，dj为支路j的管线的当量直径，Rej为支路j中介质流动时的流态判断变量雷诺数,λj为支路j的管线的沿程
阻力系数,ρbj为支路j中介质的密度；

[0020] (3)从调门曲线选取三个调门中间开度θ1,θ2,θ3，根据调门曲线选取的与开度θ1,θ2,θ3相对应的调门通流比例χ1,χ2,χ3，计算支路j的通流比例χ和实际通流系数Cvj, j＝1,
2…n，n为综合能源网中热力管网中的支路总数:

[0021] 对支路j进行判断，若支路j不带有调节阀门，则将支路j的通流比例χ置为1，若支路j带有调节阀门，则采用四次多项式拟合的方法，

[0022] 其中，aa0,aa1,aa2,aa3,aa4分别为拟合多项式系数，

[0023] 对阀门通流比例χ进行多项式拟合计算：

[0024] χ＝aa0+aa1θ+aa2θ2+aa3θ3+aa4θ4

[0025] Cvj＝χKvj

[0026] 上式中，Kvj为步骤(2)中支路j的最大通流系数，χ的取值为0‑1，θ为调节阀门的开度，θ的取值范围为0‑1；

[0027] (4)根据与支路j相对应的管线内的流体质量流量wj，利用下式，计算支路j的流导Rj：

[0028]

[0029] 上式中，ρbj为步骤(1)的支路j中流动介质的密度；

[0030] (5)根据支路j(j＝1,2…n)中流体质量流量wj、与支路j相连的上节点的海拔高度hcj、与支路j相连的下节点的海拔高度hdj，利用下式，计算与热力管网水力迭代计算过程中
支路j相对应的流量补偿项Cj：

[0031]

[0032] 式中，g为重力加速度，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；

[0033] (6)根据步骤(4)的支路j(j＝1,2…n)流导Rj，利用下式，计算得到一个用于热力管网水力计算的水力求解系数矩阵A的对角元素aii以及与节点i相对应的非对角元素 aik,
aki：

[0034]

[0035] aik＝aki＝Rjio

[0036] 式中，二维数组AI[i][ai]的行号i为节点编号，第i行的行元素为与流入编号为i的节点的支路相连的相邻节点编号，ai为与流入编号为i的节点的支路相邻的节点总数，二
维数组BI[i][ai]的行号i为节点编号，第i行的行元素为与流出编号为i的节点的支路相连
的相邻节点编号，ib为与流出编号为i的节点的支路相邻的节点总数，下标CI[i][ai+ib]为
一个二维数组，二维数组中的行号i为节点编号，第i行的行元素依次为ai个流入编号为i 
的节点的支路编号和ib个流出编号为i的节点的支路编号，iki为遍历流入编号为i的支路
的计数变量，iki＝1,2…ai，iko为遍历流出支路i的计数变量，iko＝1,2…ib，jio(jio＝CI
[i][1],CI[i][2]…CI[i][ai+ib])为与节点i的相连的支路编号，k为与节点i相连的编号
为jio的支路相连的相邻节点编号，k＝AI[i][1],…AI[i][ai],BI[i][1],…BI[i][ib]；

[0037] (7)根据步骤(5)求得的热力管网各支路的流量补偿项Cj，利用下式，计算得到一个流量补偿向量B，流量补偿向量B中的元素为bi(i＝1,2…m)：

[0038]

[0039] 式中，iki为遍历流入编号为i的节点的支路的计数变量，iki＝1,2…ai，iko为遍历流出节点i的支路的计数变量，iko＝1,2…ib；

[0040] (8)根据步骤(6)和步骤(7)计算的水力求解系数矩阵A和流量补偿向量为B的元素值，得到一个热力管网水力计算的线性方程组，利用高斯‑赛德尔迭代的方法求解该线性方
程组，得到t至t+ΔT计算时段内热力管网中各节点压力 i＝1,2…m：

[0041]

[0042] 式中， 为t至t+ΔT计算时段内热力管网中节点i的压力；

[0043] (9)记录上一计算时步中热力管网内各节点压力Pi，利用下式计算在t至t+ΔT计算时段内各节点压力与t‑ΔT至t计算时段内各节点压力差ΔPi，选取所有节点压力误差值
中的最大值作为收敛判据ε，并将更新的各节点压力 赋值给Pi：

[0044]

[0045] ε＝max{ΔP1,ΔP2…ΔPm}

[0046]

[0047] (10)根据步骤(9)的热力管网各节点的压力Pi、步骤(4)的热力管网各支路的流导Rj及步骤(5)的各支路流量补偿项Cj(j＝1,2…n)，利用下式更新支路j中的介质质量流量
wj；

[0048] wj＝Rj(Pcj‑Pdj)‑Cj

[0049] 式中，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；

[0050] (11)设置一个热力管网的水力迭代计算阈值δ，对步骤(9)的进行判断：若|ε|≥δ则返回步骤(3)，若|ε|＜δ，则进入步骤(12)；

[0051] (12)根据步骤(10)的介质质量流量wj(j＝1,2…n)及步骤(1)中求得的各支路中流动介质的比热cj,得到一个用于热力管网的热力计算的热力求解系数矩阵D，热力求解系
数矩阵D中与节点i(i＝1,2…m)相对应的非对角元素dif为：

[0052] dif＝wjocjoΔt

[0053] 式中，jo为流出编号为i的节点的支路编号,jo＝CI[i][ai+1]),CI[i][ai+2]…CI[i][ai+ib]， f为与流出节点i的支路jo相连的节点编号，f＝BI[i][1],…BI[i][ib]；

[0054] (13)根据步骤(2)的各支路对应管线的内容积VBj(j＝1,2…n)、步骤(1)的各节点处流动流体的密度ρi和比热ci，利用步骤(10)的热力管网各支路流动介质质量流量wj，利用
下式得到热力管网中各节点的当量容积Vi以及步骤(12)的用于热力管网的热力求解系数
矩阵D的对角线元素dii,i＝1,2…m:

[0055]

[0056]

[0057] 式中，Vi为编号为i的节点的当量内容积，ikio为遍历与编号为i的节点相连的支路的计数变量，ikio＝1,2…ai+ib，iki为遍历流入编号为i的节点的支路的计数变量， iki
＝1,2…ai：

[0058] (14)根据热力管网中各支路对外界的热功转换功率qj，j＝1,2…n，以及步骤(13) 的热力管网中各节点的当量容积Vi(i＝1,2…m)、步骤(1)的热力管网各节点处流动介质的
密度ρi和比热ci，利用下式得到热力管网各节点处注入的热功率qni：

[0059]

[0060] 并得到一个热力补偿向量为E，热力补偿向量为E中的元素ei：

[0061] ei＝‑Viρici+Δtqnj

[0062] 式中，iki为遍历流入编号为i的节点的支路的计数变量；

[0063] (15)根据步骤(12)‑(14)的热力计算系数矩阵D和热力补偿向量E的元素，得到用于热力管网热力计算的线性方程组，利用高斯‑赛德尔迭代方式求解该线性方程组，得到t
至t+Δt计算时段内热力网络各节点处输送介质的温度 i＝1,2…m：

[0064]

[0065] (16)根据步骤(15)的各节点流动介质的温度 i＝1,2…m，利用下式计算出热力管网热力计算的t至t+Δt计算时段内各节点中的介质温度与t‑Δt至t计算时段内各节
点介质温度差ΔTi，选取所有节点介质温度误差值中的最大值作为收敛判据ε'，并将更新
的各节点温度 赋值给Ti：

[0066]

[0067]

[0068] ε'＝max{ΔT1,ΔT2…ΔTm}；

[0069] (17)设置一个热力管网的节点温度的迭代计算阈值δ'，对步骤(16)ε'，进行判断：若|ε'|≥δ'，则返回步骤(2)，若|ε'|＜δ'，则结束计算，并将步骤(9)得到的热力管网内各
节点处介质的压力Pi、步骤(10)得到的热力管网内各支路中的介质质量流量wj、步骤 (16)
得到的各节点出流动介质的温度Ti作为综合能源网中热力管网的运行参数，实现能源网中
热力管网运行参数的动态计算。

[0070] 本发明提出的用于综合能源网中热力管网运行参数的动态计算方法，其特点和优点是：

[0071] 本发明的热力管网运行参数的动态计算方法，充分考虑了热力管网中工质物性变化、热力管网动态过程中水力工况和热力工况间的相互影响，不仅适用于热水等热力工况
对水力工况影响较小的动态过程，也能适用于导热油、蒸汽等热力工况对水力工况影响较
大的动态过程。同时，本发明方法从数学描述和表达形式上完成了统一，可实现解算方法的
统一，便于计算机编程实现，利用矩阵进行线性求解，具有很好的应用前景。

附图说明

[0072] 图1是本发明方法的流程框图。

[0073] 具体实施方法

[0074] 本发明提出的用于综合能源网中热力管网运行参数的动态计算方法，其流程框图如图 1所示，包括以下步骤：

[0075] 设定：综合能源网中热力管网中的节点总数为m，支路总数为n，热力计算步长为Δt，水力计算步长为ΔT，水力求解系数矩阵为A、流量补偿向量为B，热力求解系数矩阵为D、
热力补偿向量为E，初始化系数矩阵A、D的元素及补偿向量B、E元素为0，初始化各节点压力
Pi(i＝1,2…m)为0.1、节点温度Ti为20，根据设计工况初始化各支路流量 wj(j＝1,2…n)、
各支路与外界的热交换功率qj及各节点的海拔高度hi；根据综合能源网中热力管网的拓扑
结构,得到二维数组AI[i][ai]和BI[i][ai],其中二维数组AI[i][ai]的行号i为节点编号，
第i行的行元素为与流入节点i的支路相连的相邻节点编号，ai为与流入编号为i的节点的
支路相邻的节点总数，二维数组BI[i][ai]的行号i为节点编号，第i行的行元素为与流出编
号为i的节点的支路相连的相邻节点编号，ib为与流出编号为i的节点的支路相邻的节点总
数；

[0076] (1)根据热力管网中任意节点i内输送介质的温度Ti和介质压力Pi，i＝1,2…m，m为综合能源网中热力管网中的节点总数，利用介质不同温度对粘度的映射函数f1，根据温度、
压力求介质密度的函数f2和求介质比热的函数f3，分别计算支路j中输送介质的动力粘度
νj、密度ρbj以及第i号节点中的输送介质的密度ρi和比热ci；

[0077] νj＝f1(Tcj)

[0078] ρbj＝f2(Pcj,Tcj)

[0079] ρi＝f2(Pi,Ti)

[0080] ci＝f3(Pi,Ti)

[0081] cj＝f3(Pcj,Tcj)

[0082] 式中，j＝1,2…n，n为综合能源网中热力管网中的支路总数，cj为与支路j相连的上节点编号；

[0083] (2)根据热力管网中任意支路j内输送介质的温度Tj，j＝1,2…n，n为综合能源网中热力管网中的支路总数，利用介质不同温度对粘度的映射函数f，计算得到热力管网中任
意支路j的流体动力黏度νj，利用下式求出支路j的最大通流系数Kvj及支路j的内容积VBj；

[0084]

[0085]

[0086]

[0087]

[0088] 式中，lj为支路j的管线长度，uj为支路j的管线中输送介质的流速，dj为支路j的管线的当量直径，Rej为支路j中介质流动时的流态判断变量雷诺数,λj为支路j的管线的沿程
阻力系数,ρbj为支路j中介质的密度；

[0089] (3)从调门曲线(由阀门出厂厂家提供)选取三个调门中间开度θ1,θ2,θ3，根据调门曲线选取的与开度θ1,θ2,θ3相对应的调门通流比例χ1,χ2,χ3，计算支路j的通流比例χ和实际
通流系数Cvj：

[0090] 对支路j进行判断，若支路j不带有调节阀门，则将支路j的通流比例χ置为1，若支路j带有调节阀门，则采用四次多项式拟合的方法，

[0091] 其中，aa0,aa1,aa2,aa3,aa4分别为拟合多项式系数，

[0092] 对阀门通流比例χ进行多项式拟合计算：

[0093] χ＝aa0+aa1θ+aa2θ2+aa3θ3+aa4θ4

[0094] Cvj＝χKvj

[0095] 上式中，Kvj为步骤(2)中支路j的最大通流系数，χ的取值为0‑1，θ为调节阀门的开度，θ的取值范围为0‑1；

[0096] (4)根据与支路j相对应的管线内的流体质量流量wj，利用下式，计算支路j的流导 Rj：

[0097]

[0098] 上式中，ρbj为步骤(1)的支路j中流动介质的密度；

[0099] (5)根据支路j中流体质量流量wj、与支路j相连的上节点的海拔高度hcj、与支路j 相连的下节点的海拔高度hdj，利用下式，计算与热力管网水力迭代计算过程中支路j相对应
的流量补偿项Cj：

[0100]

[0101] 式中，g为重力加速度，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；

[0102] (6)根据步骤(4)的支路j流导Rj，利用下式，计算得到一个用于热力管网水力计算的水力求解系数矩阵A的对角元素aii以及与节点i相对应的非对角元素aik,aki：

[0103]

[0104] aik＝aki＝Rjio

[0105] 式中，二维数组AI[i][ai]的行号i为节点编号，第i行的行元素为与流入编号为i的节点的支路相连的相邻节点编号，ai为与流入编号为i的节点的支路相邻的节点总数，二
维数组BI[i][ai]的行号i为节点编号，第i行的行元素为与流出编号为i的节点的支路相连
的相邻节点编号，ib为与流出编号为i的节点的支路相邻的节点总数，下标CI[i][ai+ib]为
一个二维数组，二维数组中的行号i为节点编号，第i行的行元素依次为ai个流入编号为i 
的节点的支路编号和ib个流出编号为i的节点的支路编号，iki为遍历流入编号为i的支路
的计数变量，iki＝1,2…ai，iko为遍历流出支路i的计数变量，iko＝1,2…ib， jio(jio＝
CI[i][1],CI[i][2]…CI[i][ai+ib])为与节点i的相连的支路编号，k为与节点i相连的编
号为jio的支路相连的相邻节点编号，k＝AI[i][1],…AI[i][ai],BI[i][1],…BI[i][ib]；

[0106] (7)根据步骤(5)求得的热力管网各支路的流量补偿项Cj，利用下式，计算得到一个流量补偿向量B，流量补偿向量B中的元素为bi(i＝1,2…m)：

[0107]

[0108] 式中，iki为遍历流入编号为i的节点的支路的计数变量，iki＝1,2…ai，iko为遍历流出节点i的支路的计数变量，iko＝1,2…ib；

[0109] (8)根据步骤(6)和步骤(7)计算的水力求解系数矩阵A和流量补偿向量为B的元素值，得到一个热力管网水力计算的线性方程组，利用高斯‑赛德尔迭代的方法求解该线性方
程组，得到t至t+ΔT计算时段内热力管网中各节点压力 i＝1,2…m：

[0110]

[0111] 式中， 为t至t+ΔT计算时段内热力管网中节点i的压力；

[0112] (9)记录上一计算时步中热力管网内各节点压力Pi，利用下式计算在t至t+ΔT计算时段内各节点压力与t‑ΔT至t计算时段内各节点压力差ΔPi，选取所有节点压力误差值
中的最大值作为收敛判据ε，并将更新的各节点压力 赋值给Pi：

[0113]

[0114] ε＝max{ΔP1,ΔP2…ΔPm}

[0115]

[0116] (10)根据步骤(9)的热力管网各节点的压力Pi、步骤(4)的热力管网各支路的流导Rj及步骤(5)的各支路流量补偿项Cj，利用下式更新支路j中的介质质量流量wj；

[0117] wj＝Rj(Pcj‑Pdj)‑Cj

[0118] 式中，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；

[0119] (11)设置一个热力管网的水力迭代计算阈值δ，对步骤(9)的进行判断：若|ε|≥δ则返回步骤(3)，若|ε|＜δ，则进入步骤(12)；

[0120] (12)根据步骤(10)的介质质量流量wj(j＝1,2…n)及步骤(1)中求得的各支路中流动介质的比热cj,得到一个用于热力管网的热力计算的热力求解系数矩阵D，热力求解系
数矩阵D中与节点i相对应的非对角元素dif为：

[0121] dif＝wjocjoΔt

[0122] 式中，jo为流出编号为i的节点的支路编号,jo＝CI[i][ai+1]),CI[i][ai+2]…CI[i][ai+ib]， f为与流出节点i的支路jo相连的节点编号，f＝BI[i][1],…BI[i][ib]；

[0123] (13)根据步骤(2)的各支路对应管线的内容积VBj(j＝1,2…n)、步骤(1)的各节点处流动流体的密度ρi和比热ci，利用步骤(10)的热力管网各支路流动介质质量流量wj，利用
下式得到热力管网中各节点的当量容积Vi以及步骤(12)的用于热力管网的热力求解系数
矩阵D的对角线元素dii：

[0124]

[0125]

[0126] 式中，Vi为编号为i的节点的当量内容积，ikio为遍历与编号为i的节点相连的支路的计数变量，ikio＝1,2…ai+ib，iki为遍历流入编号为i的节点的支路的计数变量， iki
＝1,2…ai：

[0127] (14)根据热力管网中各支路对外界的热功转换功率qj，j＝1,2…n，以及步骤(13) 的热力管网中各节点的当量容积Vi(i＝1,2…m)、步骤(1)的热力管网各节点处流动介质的
密度ρi和比热ci，利用下式得到热力管网各节点处注入的热功率qni：

[0128]

[0129] 并得到一个热力补偿向量为E，热力补偿向量为E中的元素ei：

[0130] ei＝‑Viρici+Δtqnj

[0131] 式中，iki为遍历流入编号为i的节点的支路的计数变量；

[0132] (15)根据步骤(12)‑(14)的热力计算系数矩阵D和热力补偿向量E的元素，得到用于热力管网热力计算的线性方程组，利用高斯‑赛德尔迭代方式求解该线性方程组，得到t
至t+Δt计算时段内热力网络各节点处输送介质的温度 i＝1,2…m：

[0133]

[0134] (16)根据步骤(15)的各节点流动介质的温度 i＝1,2…m，利用下式计算出热力管网热力计算的t至t+Δt计算时段内各节点中的介质温度与t‑Δt至t计算时段内各节
点介质温度差ΔTi，选取所有节点介质温度误差值中的最大值作为收敛判据ε'，并将更新
的各节点温度 赋值给Ti：

[0135]

[0136]

[0137] ε'＝max{ΔT1,ΔT2…ΔTm}；

[0138] (17)设置一个热力管网的节点温度的迭代计算阈值δ'，对步骤(16)ε'，进行判断：若|ε'|≥δ'，则返回步骤(2)，若|ε'|＜δ'，则结束计算，并将步骤(9)得到的热力管网内各
节点处介质的压力Pi、步骤(10)得到的热力管网内各支路中的介质质量流量wj、步骤 (16)
得到的各节点出流动介质的温度Ti作为综合能源网中热力管网的运行参数，实现能源网中
热力管网运行参数的动态计算。