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

权利要求书

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

设定：综合能源网中热力管网中的节点总数为m，支路总数为n，热力计算步长为Δt，水力计算步长为ΔT，水力求解系数矩阵为A、流量补偿向量为B，热力求解系数矩阵为D、热力补偿向量为E，初始化系数矩阵A、D的元素及补偿向量B、E元素为0，初始化各节点压力P_i，为0.1、节点温度T_i为20，i＝1,2…m，根据设计工况初始化各支路流量w_j、各支路与外界的热交换功率q_j及各节点的海拔高度h_i，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内输送介质的温度T_i和介质压力P_i，i＝1,2…m，m为综合能源网中热力管网中的节点总数，利用介质不同温度对粘度的映射函数f₁，根据温度、压力求介质密度的函数f₂和求介质比热的函数f₃，分别计算支路j中输送介质的动力粘度ν_j、密度ρb_j以及第i号节点中的输送介质的密度ρ_i和比热c_i；

ν_j＝f₁(T_cj)

ρb_j＝f₂(P_cj,T_cj)

ρ_i＝f₂(P_i,T_i)

c_i＝f₃(P_i,T_i)

c_j＝f₃(P_cj,T_cj)

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

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

式中，l_j为支路j的管线长度，u_j为支路j的管线中输送介质的流速，d_j为支路j的管线的当量直径，R_ej为支路j中介质流动时的流态判断变量雷诺数,λ_j为支路j的管线的沿程阻力系数,ρb_j为支路j中介质的密度；

(3)从调门曲线选取三个调门中间开度θ₁,θ₂,θ₃，根据调门曲线选取的与开度θ₁,θ₂,θ₃相对应的调门通流比例χ₁,χ₂,χ₃，计算支路j的通流比例χ和实际通流系数Cv_j：

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

其中，aa₀,aa₁,aa₂,aa₃,aa₄分别为拟合多项式系数，

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

χ＝aa₀+aa₁θ+aa₂θ²+aa₃θ³+aa₄θ⁴

Cv_j＝χKv_j

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

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

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

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

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

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

a_ik＝a_ki＝R_jio

式中，二维数组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)求得的热力管网各支路的流量补偿项C_j，利用下式，计算得到一个流量补偿向量B，流量补偿向量B中的元素为b_i，i＝1,2…m：

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

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

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

(9)记录上一计算时步中热力管网内各节点压力P_i，i＝1,2…m，利用下式计算在t至t+ΔT计算时段内各节点压力与t-ΔT至t计算时段内各节点压力差ΔP_i，选取所有节点压力误差值中的最大值作为收敛判据ε，并将更新的各节点压力P_i^t+ΔT赋值给P_i：

ΔP_i＝P_i^t+ΔT-P_i

ε＝max{ΔP₁,ΔP₂…ΔP_m}

P_i＝P_i^t+ΔT

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

w_j＝R_j(P_cj-P_dj)-C_j

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

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

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

d_if＝w_joc_joΔ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)的各支路对应管线的内容积VB_j、步骤(1)的各节点处流动流体的密度ρ_i和比热c_i，i＝1,2…m，利用步骤(10)的热力管网各支路流动介质质量流量w_j，利用下式得到热力管网中各节点的当量容积V_i以及步骤(12)的用于热力管网的热力求解系数矩阵D的对角线元素d_ii：

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

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

并得到一个热力补偿向量为E，热力补偿向量为E中的元素e_i：

e_i＝-V_iρ_ic_i+Δtqn_j

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

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

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

ΔT_i＝T_i^t+Δt-T_i

T_i＝T_i^t+Δt

ε'＝max{ΔT₁,ΔT₂…ΔT_m}；

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