基于温度插值求解的管网仿真研究

2023-05-28秦书松栾积毅张文丰高建民

节能技术 2023年2期

秦书松,栾积毅,张文丰,谢敏,杜谦,高建民

(1.佳木斯大学机械工程学院,黑龙江佳木斯 154007;2.哈电发电设备国家工程研究中心,黑龙江哈尔滨 150028;3.哈尔滨工业大学能源科学与工程学院,黑龙江哈尔滨 150001)

0 引言

近年来,我国加快推动能源生产、运输等各环节的数字化智能化升级,推动能源行业的低碳转型。蒸汽供热作为能源领域的重要子行业,随着管网拓扑结构复杂性的增加、末端用户品质需求的增长和系统节能减排压力的增强,原有供热设计方法的适应性也逐渐降低[1-3]。为精准构建模型,通常要结合传热学、流体力学和建筑暖通等方面知识,这也决定了模型构建的复杂性。

应用Flowmaster、Trnsys等通用模拟软件,可快速实现管网系统模拟,辅助管网设计。但上述软件多用于石油化工等领域,且需要根据软件操作相应调整模拟过程。对蒸汽管网水力计算而言,模型建立及解法较为成熟[4]。在水力热力耦合建模方面,孙玉宝[5]通过反复迭代校正管道流量、管道平均密度和定压比热容的方式得到系统参数,经实例验证计算误差在5%以内。Wang H等[6]基于四阶Runge-Kutta算法建立蒸汽管网水力热力耦合模型,但该方法的计算速度相对较慢,不适用于结构复杂的管网。王安然等[7]在求解过程中采用隐式欧拉算法、稀疏矩阵算法求解模型非线性微分方程组,并结合Fluent和工程实际验证。陈鸿伟等[8]结合蒸汽在管道中的流动机理,提出基于蒸汽水力热力耦合计算和误差修正模型的混合建模方法,减小机理模型与实际蒸汽流动工况的误差。陈彬彬等[9]以供热网络为对象推导了统一能路理论中的水力与热力模型,并以此刻画了供热网络的支路特性和拓扑约束。算例结果表明了潮流计算方法分析供热网络的有效性。

上述研究建立在“模型-仿真计算-验证”的基础上,关键点在模型的合理性与准确性上。然而针对水力模型忽略热力工况,水力热力耦合模型线性化求解复杂的问题,本文建立温度插值策略求解的管网模型,

即在计算过程中水力模型的求解基于热力模型求解反馈的温度参数值。

1 数学模型

1.1 质能平衡

流体网络和电网都遵循质量和能量守恒定律,因此根据图论原理简化管网拓扑结构,结合“点线”、“线环”组网的方法,利用电网中的基尔霍夫相关定律可以求解蒸汽管网的运行参数[10]。

1.1.1 节点质量平衡

基尔霍夫电流定律表明,对于网络中任意节点,流入与流出该节点的代数和为零。如图1所示,可以理解为流入与流出某一节点的流量代数和为零。对于已知的有N个节点与M条管道的管网有向图G,利用基尔霍夫电流定律建立节点流量平衡方程。

图1 节点质量平衡图

(1)

式中A——流出或流入对应关联矩阵;

aij——A中的元素;

G=[G1,G2……,GM]T——管道的流量;

Gj——第j根管道的流量;

Q=[Q1,Q2,Q3,…,QN]T——节点的流量。

1.1.2 回路能量平衡

由于集中供热管网为闭合的恒定流系统,因此根据能量守恒定律,沿任一回路方向各管道压降的代数和为零。供热管网中的水力压降也同样可以类比电学中的电力压降,所以沿任意回路方向各个管道压降的代数和为零。式(2)为M阶列向量ΔH,表示系统中各管道的压降。因此结合基尔霍夫电压定律,可以得到通用的回路水力能量方程(3)

ΔH=[Δh1,Δh2,…,ΔhM]T

(2)

(3)

同理,对于管网中任一回路也均应满足温降之和为零,因此可以构建回路热力能量方程(5)。式(4)为M阶列向量ΔT,表示各个管道的温降;bsj表示管道与回路的关联信息,为基本回路矩阵B中的元素

ΔT=(Δt1,Δt2,…,ΔtM)T

(4)

(5)

1.2 热力模型

1.2.1 热力特性

蒸汽管道可近似看作单层圆筒壁,而实际的蒸汽管道是由n层不同材料组成的,因此多层圆筒壁的总热阻等于各层热阻之和。蒸汽管道的总传热系数主要与蒸汽和管道内壁的对流换热、保温层的导热、以及空气(架空)或土壤(直埋)与管道外护层的换热等因素有关。由于钢管管壁导热热阻远小于保温层导热热阻,一般可忽略不计,所以传热系数[11]可按式(6)计算,式中参数如图2中定义所示。其中保温层的导热主要与保温材质和保温层厚度有关系,在此不做详述,其他可由式(6)-式(10)推导分析。

图2 保温管道结构示意图

(6)

管道内蒸汽在热源压力影响下向前受迫流动,对流换热的实验关联式采用式(7)迪图斯-贝尔特公式计算。根据努塞尔数的定义,蒸汽与管道内壁的对流换热系数可由式(8)计算

(7)

(8)

式中Nu——努塞尔数;

Re——雷诺数;

Pr——普朗特数;

hm——蒸汽与管道内壁的对流换热系数/W·(m2·K)-1。

环境与架空管道外护层表面的对流换热系数,一般可根据具体风速按式(9)校核。土壤与地埋管道的换热系数,按福尔赫盖伊默推导的传热学理论公式计算,如式(10)所示

(9)

(10)

式中 h0——管道外护层对流换热系数/w·(m2·K)-1;

λt——土壤的导热系数/W/(m·K)-1;

w——风速/m·s-1;

h——折算地埋深度/m。

1.2.2 热力模型与求解

蒸汽各节点温度主要取决于管道和附件的沿途温降,而沿途温降主要与沿途管道的散热、蒸汽流量、管道长度以及蒸汽的定压比热有关。由此根据热平衡原理,结合节点流量平衡方程、回路能量方程和温降计算公式,建立热力求解方程组(11)[12]。如图3所示,温度插值策略求解过程中管网热力计算时的流量和压力以水力计算的结果为基础,对节点温度反复修正。不仅可以减小复杂水力条件下热力参数对水力参数的影响,提高稳定性,而且采用两个小方程组代替一个大方程组,可以加快求解速度。

图3 温度插值策略求解流程图

(11)

式中A0——由aij组成的(N-1)×M降阶关联矩阵;

F——由管道温降系数组成的M阶节点对角矩阵;

AT——矩阵A的转置矩阵;

T——节点温度向量,T=(t1,t2,t3,…,tN-1)T。

方程组(11)共有2M+N-1个方程,未知各管道的流量、温降以及各节点的温度共2M+N-1个,即方程个数与未知参数数目相等,对该方程组求解可得节点温度方程(13)

[A·(F|G|)-1·AT]T+Q=0

(12)

T=(A·(F|G|)-1·AT)-1Q

(13)

通过上述热力计算模型求解分析可知,在节点温度及管道温降求解过程中管道流量对计算结果有着较为重要的影响,为使系统计算得到的结果更加符合实际并减小误差,应首先对系统流量进行分配。

1.3 水力模型

1.3.1 水力特性

供热管道的压力损失等于始、末两断面间沿程和局部阻力损失之和,如式(14)所示。在单位长度范围内,假设管道摩擦阻力系数和管内蒸汽密度不随管道长度变化,沿程阻力损失可由Darcy公式(15)计算

ΔHh=ΔHl+ΔHf

(14)

(15)

式中 ΔHh——压力损失/Pa;

ΔHl——沿程阻力损失/Pa;

ΔHf——局部阻力损失/Pa;

L——长度/m;

λ——阻力系数;

ρ——密度/kg·m-3;

V——速度/m·s-1。

大部分蒸汽的流动都属于紊流过渡区,如式(16)所示,Colebrook-White公式给出了管道当量粗糙度、雷诺数、阻力系数的关系,适用于工业管道的三个阻力区,又称为紊流λ的综合公式。阻力系数代入Darcy公式,可以得到目前理论上最精确的压力损失计算结果。蒸汽流经弯头、三通等部件所产生的局部阻力损失,通常采用局部阻力系数和当量长度近似方法计算。结合管道内蒸汽流速计算公式(17),可以得到局部阻力损失计算公式(18)

(16)

(17)

(18)

式中ε——管道相对粗糙度/m;

Le——当量长度/m。

假定管道内的凝结水均在节点处生成排出,蒸汽流量不变,管道内流体的压降和流量的关系服从二次幂规律。蒸汽管网稳态水力计算压力损失数学表达式可以归纳为式(19)

(19)

式中S——管道阻力特性系数/Pa·(kg·h)-2。

1.3.2 水力模型与求解

对于已知的管网系统,参照电路理论中的基尔霍夫定律,结合阻力损失公式、节点流量平衡和回路能量平衡方程等建立通用的蒸汽管网流量分配方程组(20)[13]

(20)

式中S=diag(s1,s2,s3,…,sM)——管道阻力特性系数对角矩阵;

H=(H1,H2,H3,…HN-1)T——节点压力向量;

|G|=diag(|G1|,|G2|,|G3|,…|GM|)——管道流量的对角矩阵。

由方程组(20)可知该管道流量分配模型中管道流量满足节点流量平衡方程,同时由方程组(11)可知热力计算过程中管道流量同样满足节点流量平衡方程。因此,热力计算时的管道流量采用该模型计算得出的结果的方法可行。式(21)和(22)为方程组(21)化简求解过程

[A·(S|G|)-1·AT]H+Q=0

(21)

H=(A·(S|G|)-1·AT)-1Q

(22)

2 模型线性化求解与案例分析

2.1 模型线性化求解

上述模型仅包含流量、压力和温度三个基础变量,式(12)和(21)中(A·(F|G|)-1·AT)-1、(A·(S|G|)-1·AT)-1实际为N·N阶管网系数矩阵。根据节点分类的不同,划分管网系数矩阵,可以得到控制性方程组(23),目前被认为是最简洁、高效的管网计算模型[14]

(23)

式中A11、0、A21、A22——系数矩阵划分后的分块矩阵;

A10=[M;N-n0]——管道与已知压力节点的连接矩阵;

A21——阻力特性系数对角矩阵S;

QT=[q1,q2,…,qN-n0]——已知节点需求量;

HT=[H1,H2,…,Hn0]——未知节点压力。

为方便求解,将非线性方程进行线性化处理,用求解线性方程组的方法逐步逼近非线性方程组的解,并利用迭代的方法消除误差。具体数值求解过程为式(24)～式(26)。利用Newton-Jacobi迭代方法可求取关于管道流量与节点压力的修正值,基于此方程,系统模型可由式(26)逐步迭代求解

(24)

(25)

(26)

以上为管网系统水力计算的主要过程,将管道压力降与管道流量之间的非线性关系转化为线性关系,迭代修正参数与基值之间的误差,使前后两次的值不断接近,获得系统的压力和流量分布。针对热力计算,可以采用相似的数值求解方法。

2.2 案例分析

为验证本文所建立模型的有效性,对某工业园区进行仿真分析。根据供热系统结构简化拓扑关系,如图4所示:由29个节点、28条管道组成,每条管道末端为热用户,主要有16家用热企业。该蒸汽管网为枝状,全长约27 km,管径范围从DN70到 DN600。热源厂供给压力为1.6 MPa,温度为250 ℃。管网具体设计参数如表1所示。

图4 园区供热管网规划示意图

图5 节点运行和仿真数据对比

图6 压力和温度误差分析

表1 园区管网设计参数

由于该工业园区尚未对凝结水量进行监测,提供流量参数仅有热源供给量和热用户需求量,所建立的一维稳态模型假设蒸汽过热输送且为单一均相,无法对管道凝结水量进行评估。因此不再对管道蒸汽和凝结水流量讨论,仿真结果如表2所示。

表2 管网系统仿真结果

续表2

通过以上图表分析可知,在大部分节点处,仿真计算结果与实际运行数据吻合度较好,温度和压力误差大多在5%以内。节点压力和温度,在上游靠近热源点处与实际运行数据较为接近,但随着与热源点距离的增加,相对误差不断增大,节点F2、G、f2、F、f的误差相对较大。其中F点压力误差最大为 5.74%,a点压力误差最小为 0.193%,平均压力误差为 3.05%;f点温度误差最大为 5.45%,C11点温度误差最小为 1.10%,平均温度误差为2.69%。这是由于节点参数是通过与热源点的差值迭代求得,因此随着距离的增加计算误差不断积累变大。