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

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

Info

Publication number
CN111611690A
CN111611690A CN202010307723.1A CN202010307723A CN111611690A CN 111611690 A CN111611690 A CN 111611690A CN 202010307723 A CN202010307723 A CN 202010307723A CN 111611690 A CN111611690 A CN 111611690A
Authority
CN
China
Prior art keywords
node
branch
pipe network
medium
network
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010307723.1A
Other languages
English (en)
Other versions
CN111611690B (zh
Inventor
胡宪法
张树卿
唐绍普
彭振
刘宁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tsinghua University
Original Assignee
Tsinghua University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tsinghua University filed Critical Tsinghua University
Priority to CN202010307723.1A priority Critical patent/CN111611690B/zh
Publication of CN111611690A publication Critical patent/CN111611690A/zh
Application granted granted Critical
Publication of CN111611690B publication Critical patent/CN111611690B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Systems or methods specially adapted for specific business sectors, e.g. utilities or tourism
    • G06Q50/06—Electricity, gas or water supply
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation

Abstract

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

Description

一种用于综合能源网中热力管网运行参数的动态计算方法
技术领域
本发明涉及一种用于综合能源网中热力管网运行参数的动态计算方法，属于综合能源网数字仿真技术领域。
背景技术
综合能源网覆盖了热力侧子网、机电侧子网、燃料侧子网等不同类型能量网络及其具有不同响应时间的设备部件等，热力管网的动态仿真技术在综合能源系统仿真中起到越来越重要的作用。
目前的热力管网动态仿真建模主要从流体动力学的角度完成水力工况的建模，针对复杂的热力工况、工质物性变化、相变等因素缺乏统一的处理手段，且当热力管网和电力网联立解算的时候因建模方法和数学表达不统一无法统一其算法，程序实现困难，为综合能源系统仿真建模带来困扰。如导热油输送管网在导热油温度上升时，未考虑到导热油的物性影响，随着温度上升，导热油因为动力粘度的变化引起流阻的变化，热水输送过程中，未考虑因热水温度上升可能引起相变的影响，忽略了热力工况对水力工况的影响，孤立的对待热力工况和水力工况。
发明内容
本发明的目的提出一种用于综合能源网中热力管网运行参数的动态计算方法，借鉴机电侧子网机-网建模的方法，根据传统热力学理论从“广义电路”的角度对热力管网模型进行统一描述和数学表达，完成综合能源系统动态仿真模型求解的算法统一和程序实现。
本发明提出的用于综合能源网中热力管网运行参数的动态计算方法，包括以下步骤：
设定：综合能源网中热力管网中的节点总数为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的节点的支路相邻的节点总数；
(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；
Figure BDA0002456369170000021
Figure BDA0002456369170000022
Figure BDA0002456369170000031
Figure BDA0002456369170000032
式中，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＝1,2…n，n为综合能源网中热力管网中的支路总数:
对支路j进行判断，若支路j不带有调节阀门，则将支路j的通流比例χ置为1，若支路j带有调节阀门，则采用四次多项式拟合的方法，
Figure BDA0002456369170000033
其中，aa0,aa1,aa2,aa3,aa4分别为拟合多项式系数，
对阀门通流比例χ进行多项式拟合计算：
χ＝aa0+aa1θ+aa2θ2+aa3θ3+aa4θ4
Cvj＝χKvj
上式中，Kvj为步骤(2)中支路j的最大通流系数，χ的取值为0-1，θ为调节阀门的开度，θ的取值范围为0-1；
(4)根据与支路j相对应的管线内的流体质量流量wj，利用下式，计算支路j的流导Rj：
Figure BDA0002456369170000041
上式中，ρbj为步骤(1)的支路j中流动介质的密度；
(5)根据支路j(j＝1,2…n)中流体质量流量wj、与支路j相连的上节点的海拔高度hcj、与支路j相连的下节点的海拔高度hdj，利用下式，计算与热力管网水力迭代计算过程中支路j相对应的流量补偿项Cj：
Figure BDA0002456369170000042
式中，g为重力加速度，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；
(6)根据步骤(4)的支路j(j＝1,2…n)流导Rj，利用下式，计算得到一个用于热力管网水力计算的水力求解系数矩阵A的对角元素aii以及与节点i相对应的非对角元素aik,aki：
Figure BDA0002456369170000043
aik＝aki＝Rjio
式中，二维数组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]；
(7)根据步骤(5)求得的热力管网各支路的流量补偿项Cj，利用下式，计算得到一个流量补偿向量B，流量补偿向量B中的元素为bi(i＝1,2…m)：
Figure BDA0002456369170000051
式中，iki为遍历流入编号为i的节点的支路的计数变量，iki＝1,2…ai，iko为遍历流出节点i的支路的计数变量，iko＝1,2…ib；
(8)根据步骤(6)和步骤(7)计算的水力求解系数矩阵A和流量补偿向量为B的元素值，得到一个热力管网水力计算的线性方程组，利用高斯-赛德尔迭代的方法求解该线性方程组，得到t至t+ΔT计算时段内热力管网中各节点压力
Figure BDA0002456369170000052
i＝1,2…m：
Figure BDA0002456369170000053
式中，
Figure BDA0002456369170000054
为t至t+ΔT计算时段内热力管网中节点i的压力；
(9)记录上一计算时步中热力管网内各节点压力Pi，利用下式计算在t至t+ΔT计算时段内各节点压力与t-ΔT至t计算时段内各节点压力差ΔPi，选取所有节点压力误差值中的最大值作为收敛判据ε，并将更新的各节点压力
Figure BDA0002456369170000055
赋值给Pi：
Figure BDA0002456369170000056
ε＝max{ΔP1,ΔP2…ΔPm}
Figure BDA0002456369170000057
(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(j＝1,2…n)及步骤(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(j＝1,2…n)、步骤(1)的各节点处流动流体的密度ρi和比热ci，利用步骤(10)的热力管网各支路流动介质质量流量wj，利用下式得到热力管网中各节点的当量容积Vi以及步骤(12)的用于热力管网的热力求解系数矩阵D的对角线元素dii,i＝1,2…m:
Figure BDA0002456369170000061
Figure BDA0002456369170000062
式中，Vi为编号为i的节点的当量内容积，ikio为遍历与编号为i的节点相连的支路的计数变量，ikio＝1,2…ai+ib，iki为遍历流入编号为i的节点的支路的计数变量，iki＝1,2…ai：
(14)根据热力管网中各支路对外界的热功转换功率qj，j＝1,2…n，以及步骤(13)的热力管网中各节点的当量容积Vi(i＝1,2…m)、步骤(1)的热力管网各节点处流动介质的密度ρi和比热ci，利用下式得到热力管网各节点处注入的热功率qni：
Figure BDA0002456369170000071
并得到一个热力补偿向量为E，热力补偿向量为E中的元素ei：
ei＝-Viρici+Δtqnj
式中，iki为遍历流入编号为i的节点的支路的计数变量；
(15)根据步骤(12)-(14)的热力计算系数矩阵D和热力补偿向量E的元素，得到用于热力管网热力计算的线性方程组，利用高斯-赛德尔迭代方式求解该线性方程组，得到t至t+Δt计算时段内热力网络各节点处输送介质的温度
Figure BDA0002456369170000072
i＝1,2…m：
Figure BDA0002456369170000073
(16)根据步骤(15)的各节点流动介质的温度
Figure BDA0002456369170000074
i＝1,2…m，利用下式计算出热力管网热力计算的t至t+Δt计算时段内各节点中的介质温度与t-Δt至t计算时段内各节点介质温度差ΔTi，选取所有节点介质温度误差值中的最大值作为收敛判据ε'，并将更新的各节点温度
Figure BDA0002456369170000075
赋值给Ti：
Figure BDA0002456369170000076
Figure BDA0002456369170000077
ε'＝max{ΔT1,ΔT2…ΔTm}；
(17)设置一个热力管网的节点温度的迭代计算阈值δ'，对步骤(16)ε'，进行判断：若|ε'|≥δ'，则返回步骤(2)，若|ε'|＜δ'，则结束计算，并将步骤(9)得到的热力管网内各节点处介质的压力Pi、步骤(10)得到的热力管网内各支路中的介质质量流量wj、步骤(16)得到的各节点出流动介质的温度Ti作为综合能源网中热力管网的运行参数，实现能源网中热力管网运行参数的动态计算。
本发明提出的用于综合能源网中热力管网运行参数的动态计算方法，其特点和优点是：
本发明的热力管网运行参数的动态计算方法，充分考虑了热力管网中工质物性变化、热力管网动态过程中水力工况和热力工况间的相互影响，不仅适用于热水等热力工况对水力工况影响较小的动态过程，也能适用于导热油、蒸汽等热力工况对水力工况影响较大的动态过程。同时，本发明方法从数学描述和表达形式上完成了统一，可实现解算方法的统一，便于计算机编程实现，利用矩阵进行线性求解，具有很好的应用前景。
附图说明
图1是本发明方法的流程框图。
具体实施方法
本发明提出的用于综合能源网中热力管网运行参数的动态计算方法，其流程框图如图1所示，包括以下步骤：
设定：综合能源网中热力管网中的节点总数为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的节点的支路相邻的节点总数；
(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为综合能源网中热力管网中的支路总数，利用介质不同温度对粘度的映射函数f，计算得到热力管网中任意支路j的流体动力黏度νj，利用下式求出支路j的最大通流系数Kvj及支路j的内容积VBj；
Figure BDA0002456369170000091
Figure BDA0002456369170000092
Figure BDA0002456369170000093
Figure BDA0002456369170000094
式中，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带有调节阀门，则采用四次多项式拟合的方法，
Figure BDA0002456369170000101
其中，aa0,aa1,aa2,aa3,aa4分别为拟合多项式系数，
对阀门通流比例χ进行多项式拟合计算：
χ＝aa0+aa1θ+aa2θ2+aa3θ3+aa4θ4
Cvj＝χKvj
上式中，Kvj为步骤(2)中支路j的最大通流系数，χ的取值为0-1，θ为调节阀门的开度，θ的取值范围为0-1；
(4)根据与支路j相对应的管线内的流体质量流量wj，利用下式，计算支路j的流导Rj：
Figure BDA0002456369170000102
上式中，ρbj为步骤(1)的支路j中流动介质的密度；
(5)根据支路j中流体质量流量wj、与支路j相连的上节点的海拔高度hcj、与支路j相连的下节点的海拔高度hdj，利用下式，计算与热力管网水力迭代计算过程中支路j相对应的流量补偿项Cj：
Figure BDA0002456369170000103
式中，g为重力加速度，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；
(6)根据步骤(4)的支路j流导Rj，利用下式，计算得到一个用于热力管网水力计算的水力求解系数矩阵A的对角元素aii以及与节点i相对应的非对角元素aik,aki：
Figure BDA0002456369170000111
aik＝aki＝Rjio
式中，二维数组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]；
(7)根据步骤(5)求得的热力管网各支路的流量补偿项Cj，利用下式，计算得到一个流量补偿向量B，流量补偿向量B中的元素为bi(i＝1,2…m)：
Figure BDA0002456369170000112
式中，iki为遍历流入编号为i的节点的支路的计数变量，iki＝1,2…ai，iko为遍历流出节点i的支路的计数变量，iko＝1,2…ib；
(8)根据步骤(6)和步骤(7)计算的水力求解系数矩阵A和流量补偿向量为B的元素值，得到一个热力管网水力计算的线性方程组，利用高斯-赛德尔迭代的方法求解该线性方程组，得到t至t+ΔT计算时段内热力管网中各节点压力
Figure BDA0002456369170000113
i＝1,2…m：
Figure BDA0002456369170000121
式中，
Figure BDA0002456369170000122
为t至t+ΔT计算时段内热力管网中节点i的压力；
(9)记录上一计算时步中热力管网内各节点压力Pi，利用下式计算在t至t+ΔT计算时段内各节点压力与t-ΔT至t计算时段内各节点压力差ΔPi，选取所有节点压力误差值中的最大值作为收敛判据ε，并将更新的各节点压力
Figure BDA0002456369170000123
赋值给Pi：
Figure BDA0002456369170000124
ε＝max{ΔP1,ΔP2…ΔPm}
Figure BDA0002456369170000125
(10)根据步骤(9)的热力管网各节点的压力Pi、步骤(4)的热力管网各支路的流导Rj及步骤(5)的各支路流量补偿项Cj，利用下式更新支路j中的介质质量流量wj；
wj＝Rj(Pcj-Pdj)-Cj
式中，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；
(11)设置一个热力管网的水力迭代计算阈值δ，对步骤(9)的进行判断：若|ε|≥δ则返回步骤(3)，若|ε|＜δ，则进入步骤(12)；
(12)根据步骤(10)的介质质量流量wj(j＝1,2…n)及步骤(1)中求得的各支路中流动介质的比热cj,得到一个用于热力管网的热力计算的热力求解系数矩阵D，热力求解系数矩阵D中与节点i相对应的非对角元素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(j＝1,2…n)、步骤(1)的各节点处流动流体的密度ρi和比热ci，利用步骤(10)的热力管网各支路流动介质质量流量wj，利用下式得到热力管网中各节点的当量容积Vi以及步骤(12)的用于热力管网的热力求解系数矩阵D的对角线元素dii：
Figure BDA0002456369170000131
Figure BDA0002456369170000132
式中，Vi为编号为i的节点的当量内容积，ikio为遍历与编号为i的节点相连的支路的计数变量，ikio＝1,2…ai+ib，iki为遍历流入编号为i的节点的支路的计数变量，iki＝1,2…ai：
(14)根据热力管网中各支路对外界的热功转换功率qj，j＝1,2…n，以及步骤(13)的热力管网中各节点的当量容积Vi(i＝1,2…m)、步骤(1)的热力管网各节点处流动介质的密度ρi和比热ci，利用下式得到热力管网各节点处注入的热功率qni：
Figure BDA0002456369170000133
并得到一个热力补偿向量为E，热力补偿向量为E中的元素ei：
ei＝-Viρici+Δtqnj
式中，iki为遍历流入编号为i的节点的支路的计数变量；
(15)根据步骤(12)-(14)的热力计算系数矩阵D和热力补偿向量E的元素，得到用于热力管网热力计算的线性方程组，利用高斯-赛德尔迭代方式求解该线性方程组，得到t至t+Δt计算时段内热力网络各节点处输送介质的温度
Figure BDA0002456369170000134
i＝1,2…m：
Figure BDA0002456369170000135
(16)根据步骤(15)的各节点流动介质的温度
Figure BDA0002456369170000141
i＝1,2…m，利用下式计算出热力管网热力计算的t至t+Δt计算时段内各节点中的介质温度与t-Δt至t计算时段内各节点介质温度差ΔTi，选取所有节点介质温度误差值中的最大值作为收敛判据ε'，并将更新的各节点温度
Figure BDA0002456369170000142
赋值给Ti：
Figure BDA0002456369170000143
Figure BDA0002456369170000144
ε'＝max{ΔT1,ΔT2…ΔTm}；
(17)设置一个热力管网的节点温度的迭代计算阈值δ'，对步骤(16)ε'，进行判断：若|ε'|≥δ'，则返回步骤(2)，若|ε'|＜δ'，则结束计算，并将步骤(9)得到的热力管网内各节点处介质的压力Pi、步骤(10)得到的热力管网内各支路中的介质质量流量wj、步骤(16)得到的各节点出流动介质的温度Ti作为综合能源网中热力管网的运行参数，实现能源网中热力管网运行参数的动态计算。

Claims (1)

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；
Figure FDA0002456369160000021
Figure FDA0002456369160000022
Figure FDA0002456369160000023
Figure FDA0002456369160000024
式中，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带有调节阀门，则采用四次多项式拟合的方法，
Figure FDA0002456369160000025
其中，aa0,aa1,aa2,aa3,aa4分别为拟合多项式系数，
对阀门通流比例χ进行多项式拟合计算：
χ＝aa0+aa1θ+aa2θ2+aa3θ3+aa4θ4
Cvj＝χKvj
上式中，Kvj为步骤(2)中支路j的最大通流系数，χ为调门的铜留比例，取值为0-1，θ为调节阀门的实际开度，θ的取值范围为0-1；
(4)根据与支路j相对应的管线内的流体质量流量wj，利用下式，计算支路j的流导Rj：
Figure FDA0002456369160000031
上式中，ρbj为步骤(1)的支路j中流动介质的密度；
(5)根据支路j中流体质量流量wj、与支路j相连的上节点的海拔高度hcj、与支路j相连的下节点的海拔高度hdj，利用下式，计算与热力管网水力迭代计算过程中支路j相对应的流量补偿项Cj：
Figure FDA0002456369160000032
式中，g为重力加速度，cj为与支路j相连的上节点编号，dj为与支路j相连的下节点编号；
(6)根据步骤(4)的支路j流导Rj，利用下式，计算得到一个用于热力管网水力计算的水力求解系数矩阵A的对角元素aii以及与节点i相对应的非对角元素aik,aki：
Figure FDA0002456369160000033
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：
Figure FDA0002456369160000041
式中，iki为遍历流入编号为i的节点的支路的计数变量，iki＝1,2…ai，iko为遍历流出节点i的支路的计数变量，iko＝1,2…ib；
(8)根据步骤(6)和步骤(7)计算的水力求解系数矩阵A和流量补偿向量为B的元素值，得到一个热力管网水力计算的线性方程组，利用高斯-赛德尔迭代的方法求解该线性方程组，得到t至t+ΔT计算时段内热力管网中各节点压力Pi t+ΔT，i＝1,2…m：
Figure FDA0002456369160000042
式中，Pi t+ΔT为t至t+ΔT计算时段内热力管网中节点i的压力；
(9)记录上一计算时步中热力管网内各节点压力Pi，i＝1,2…m，利用下式计算在t至t+ΔT计算时段内各节点压力与t-ΔT至t计算时段内各节点压力差ΔPi，选取所有节点压力误差值中的最大值作为收敛判据ε，并将更新的各节点压力Pi t+ΔT赋值给Pi：
ΔPi＝Pi t+ΔT-Pi
ε＝max{ΔP1,ΔP2…ΔPm}
Pi＝Pi t+ΔT
(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：
Figure FDA0002456369160000051
Figure FDA0002456369160000052
式中，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：
Figure FDA0002456369160000061
并得到一个热力补偿向量为E，热力补偿向量为E中的元素ei：
ei＝-Viρici+Δtqnj
式中，iki为遍历流入编号为i的节点的支路的计数变量；
(15)根据步骤(12)-(14)的热力计算系数矩阵D和热力补偿向量E的元素，得到用于热力管网热力计算的线性方程组，利用高斯-赛德尔迭代方式求解该线性方程组，得到t至t+Δt计算时段内热力网络各节点处输送介质的温度Ti t+Δt：
Figure FDA0002456369160000062
(16)根据步骤(15)的各节点流动介质的温度Ti t+Δt，i＝1,2…m，利用下式计算出热力管网热力计算的t至t+Δt计算时段内各节点中的介质温度与t-Δt至t计算时段内各节点介质温度差ΔTi，选取所有节点介质温度误差值中的最大值作为收敛判据ε'，并将更新的各节点温度Ti t+Δt赋值给Ti：
ΔTi＝Ti t+Δt-Ti
Ti＝Ti t+Δt
ε'＝max{ΔT1,ΔT2…ΔTm}；
(17)设置一个热力管网的节点温度的迭代计算阈值δ'，对步骤(16)ε'，进行判断：若|ε'|≥δ'，则返回步骤(2)，若|ε'|＜δ'，则结束计算，并将步骤(9)得到的热力管网内各节点处介质的压力Pi、步骤(10)得到的热力管网内各支路中的介质质量流量wj、步骤(16)得到的各节点出流动介质的温度Ti作为综合能源网中热力管网的运行参数，实现能源网中热力管网运行参数的动态计算。
CN202010307723.1A 2020-04-17 2020-04-17 一种用于综合能源网中热力管网运行参数的动态计算方法 Active CN111611690B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010307723.1A CN111611690B (zh) 2020-04-17 2020-04-17 一种用于综合能源网中热力管网运行参数的动态计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010307723.1A CN111611690B (zh) 2020-04-17 2020-04-17 一种用于综合能源网中热力管网运行参数的动态计算方法

Publications (2)

Publication Number Publication Date
CN111611690A true CN111611690A (zh) 2020-09-01
CN111611690B CN111611690B (zh) 2021-06-22

Family

ID=72195951

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010307723.1A Active CN111611690B (zh) 2020-04-17 2020-04-17 一种用于综合能源网中热力管网运行参数的动态计算方法

Country Status (1)

Country Link
CN (1) CN111611690B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114893715A (zh) * 2022-04-02 2022-08-12 安徽宇航派蒙健康科技股份有限公司 加热控制方法及其装置、系统、计算机设备和存储介质

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007333246A (ja) * 2006-06-12 2007-12-27 Nissan Motor Co Ltd 可変コンダクタンスヒートパイプおよびこれを用いた温度調節装置
CN101493861A (zh) * 2009-01-15 2009-07-29 哈尔滨工业大学 一种自然循环热水锅炉水动力数值计算方法
CN106682369A (zh) * 2017-02-27 2017-05-17 常州英集动力科技有限公司 供热管网水力仿真模型辨识修正方法及系统、操作方法
CN109255466A (zh) * 2018-07-20 2019-01-22 清华大学 一种基于多工况量测的热网稳态运行参数估计方法
CN109472413A (zh) * 2018-11-14 2019-03-15 南方电网科学研究院有限责任公司 考虑热管网传输特性的园区综合能源系统优化调度方法
CN110348602A (zh) * 2019-06-06 2019-10-18 国网浙江省电力有限公司经济技术研究院 计及天然气管网和热力管网特性的综合能源系统优化方法
CN110502791A (zh) * 2019-07-22 2019-11-26 清华大学 基于能源集线器的综合能源系统稳态建模方法
CN110717226A (zh) * 2019-09-30 2020-01-21 国网浙江杭州市萧山区供电有限公司 一种考虑能源网拓扑特性的区域综合能源系统布局规划方法
CN110909468A (zh) * 2019-11-22 2020-03-24 清华大学 一种用于综合能源网动态混合仿真的热电接口交互方法
CN111008468A (zh) * 2019-11-29 2020-04-14 上海科梁信息工程股份有限公司 综合能源能量管理系统的测试方法及其测试系统

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007333246A (ja) * 2006-06-12 2007-12-27 Nissan Motor Co Ltd 可変コンダクタンスヒートパイプおよびこれを用いた温度調節装置
CN101493861A (zh) * 2009-01-15 2009-07-29 哈尔滨工业大学 一种自然循环热水锅炉水动力数值计算方法
CN106682369A (zh) * 2017-02-27 2017-05-17 常州英集动力科技有限公司 供热管网水力仿真模型辨识修正方法及系统、操作方法
CN109255466A (zh) * 2018-07-20 2019-01-22 清华大学 一种基于多工况量测的热网稳态运行参数估计方法
CN109472413A (zh) * 2018-11-14 2019-03-15 南方电网科学研究院有限责任公司 考虑热管网传输特性的园区综合能源系统优化调度方法
CN110348602A (zh) * 2019-06-06 2019-10-18 国网浙江省电力有限公司经济技术研究院 计及天然气管网和热力管网特性的综合能源系统优化方法
CN110502791A (zh) * 2019-07-22 2019-11-26 清华大学 基于能源集线器的综合能源系统稳态建模方法
CN110717226A (zh) * 2019-09-30 2020-01-21 国网浙江杭州市萧山区供电有限公司 一种考虑能源网拓扑特性的区域综合能源系统布局规划方法
CN110909468A (zh) * 2019-11-22 2020-03-24 清华大学 一种用于综合能源网动态混合仿真的热电接口交互方法
CN111008468A (zh) * 2019-11-29 2020-04-14 上海科梁信息工程股份有限公司 综合能源能量管理系统的测试方法及其测试系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
HAMID REZA MASSRUR ET AL.: "Fast Decomposed Energy Flow in Large-Scale", 《IEEE TRANSACTIONS ON SUSTAINABLE ENERGY》 *
张树卿等: "多能微网阶段化多模式混合仿真关键技术综述", 《计算机研究与发展》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114893715A (zh) * 2022-04-02 2022-08-12 安徽宇航派蒙健康科技股份有限公司 加热控制方法及其装置、系统、计算机设备和存储介质
CN114893715B (zh) * 2022-04-02 2023-11-21 安徽宇航派蒙健康科技股份有限公司 加热控制方法及其装置、系统、计算机设备和存储介质

Also Published As

Publication number Publication date
CN111611690B (zh) 2021-06-22

Similar Documents

Publication Publication Date Title
CN105910169B (zh) 基于机理模型预测控制的城市供热系统热网调节方法及系统
CN107563674B (zh) 一种考虑管道动态特性的电-热耦合系统状态估计方法
Islamoglu et al. Heat transfer analysis using ANNs with experimental data for air flowing in corrugated channels
CN108920866B (zh) 基于滚动时域估计理论的热网动态调节运行参数估计方法
CN112417662B (zh) 一种实现集中供热管网系统动态水力优化的方法
WO2016156556A1 (en) Method and system for determining characteristic parameters of a hydraulic network
WO2015171196A1 (en) Virtual flow measurement system
CN111611690B (zh) 一种用于综合能源网中热力管网运行参数的动态计算方法
CN110135013B (zh) 一种构建供热管网动态热力模型的整合方法
CN111898816A (zh) 一种综合能源系统动态状态估计方法
CN114202111A (zh) 基于粒子群优化bp神经网络的电子膨胀阀流量特性预测
CN113111515B (zh) 一种综合能源系统的统一建模方法
CN114169239A (zh) 基于遗传算法的蒸汽管网阻力系数和换热系数辨识方法
Li et al. Study of distributed model predictive control for radiant floor heating systems with hydraulic coupling
CN105447256B (zh) 一种增强激励仿真遗传优化方法
Ermis ANN modeling of compact heat exchangers
CN111310343B (zh) 一种用于综合能源系统调度的供热网络热路建模方法
Hanafi et al. Heat exchanger's shell and tube modeling for intelligent control design
CN111783309A (zh) 基于内部守恒的蒸汽供热网络动态仿真方法
CN116611706A (zh) 基于多能源主体的动态碳排放因子测算方法
CN114117819B (zh) 一种热蒸汽网络稳态仿真方法
CN115062555A (zh) 一种基于非平衡节点㶲的综合能源系统㶲流直接计算方法
CN113111555B (zh) 一种基于叠加解耦法的质调节热力系统能流快速计算方法
CN114239199A (zh) 一种考虑凝结水的蒸汽管网动态仿真方法
CN111274668A (zh) 一种提升集中供热管网时序数据质量方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant