秦書松,欒積毅,張文豐,謝 敏,杜 謙,高建民

(1.佳木斯大學 機械工程學院,黑龍江 佳木斯 154007;2.哈電發(fā)電設備國家工程研究中心,黑龍江 哈爾濱 150028;3.哈爾濱工業(yè)大學 能源科學與工程學院,黑龍江 哈爾濱 150001)

0 引言

近年來,我國加快推動能源生產、運輸?shù)雀鳝h(huán)節(jié)的數(shù)字化智能化升級,推動能源行業(yè)的低碳轉型。蒸汽供熱作為能源領域的重要子行業(yè),隨著管網(wǎng)拓撲結構復雜性的增加、末端用戶品質需求的增長和系統(tǒng)節(jié)能減排壓力的增強,原有供熱設計方法的適應性也逐漸降低[1-3]。為精準構建模型,通常要結合傳熱學、流體力學和建筑暖通等方面知識,這也決定了模型構建的復雜性。

應用Flowmaster、Trnsys等通用模擬軟件,可快速實現(xiàn)管網(wǎng)系統(tǒng)模擬,輔助管網(wǎng)設計。但上述軟件多用于石油化工等領域,且需要根據(jù)軟件操作相應調整模擬過程。對蒸汽管網(wǎng)水力計算而言,模型建立及解法較為成熟[4]。在水力熱力耦合建模方面,孫玉寶[5]通過反復迭代校正管道流量、管道平均密度和定壓比熱容的方式得到系統(tǒng)參數(shù),經實例驗證計算誤差在5%以內。Wang H等[6]基于四階Runge-Kutta算法建立蒸汽管網(wǎng)水力熱力耦合模型,但該方法的計算速度相對較慢,不適用于結構復雜的管網(wǎng)。王安然等[7]在求解過程中采用隱式歐拉算法、稀疏矩陣算法求解模型非線性微分方程組,并結合Fluent和工程實際驗證。陳鴻偉等[8]結合蒸汽在管道中的流動機理,提出基于蒸汽水力熱力耦合計算和誤差修正模型的混合建模方法,減小機理模型與實際蒸汽流動工況的誤差。陳彬彬等[9]以供熱網(wǎng)絡為對象推導了統(tǒng)一能路理論中的水力與熱力模型,并以此刻畫了供熱網(wǎng)絡的支路特性和拓撲約束。算例結果表明了潮流計算方法分析供熱網(wǎng)絡的有效性。

上述研究建立在“模型-仿真計算-驗證”的基礎上,關鍵點在模型的合理性與準確性上。然而針對水力模型忽略熱力工況,水力熱力耦合模型線性化求解復雜的問題,本文建立溫度插值策略求解的管網(wǎng)模型,

即在計算過程中水力模型的求解基于熱力模型求解反饋的溫度參數(shù)值。

1 數(shù)學模型

1.1 質能平衡

流體網(wǎng)絡和電網(wǎng)都遵循質量和能量守恒定律,因此根據(jù)圖論原理簡化管網(wǎng)拓撲結構,結合“點線”、“線環(huán)”組網(wǎng)的方法,利用電網(wǎng)中的基爾霍夫相關定律可以求解蒸汽管網(wǎng)的運行參數(shù)[10]。

1.1.1 節(jié)點質量平衡

基爾霍夫電流定律表明,對于網(wǎng)絡中任意節(jié)點,流入與流出該節(jié)點的代數(shù)和為零。如圖1所示,可以理解為流入與流出某一節(jié)點的流量代數(shù)和為零。對于已知的有N個節(jié)點與M條管道的管網(wǎng)有向圖G,利用基爾霍夫電流定律建立節(jié)點流量平衡方程。

圖1 節(jié)點質量平衡圖

(1)

式中A——流出或流入對應關聯(lián)矩陣;

aij——A中的元素;

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

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

Q=[Q1,Q2,Q3,…,QN]T——節(jié)點的流量。

1.1.2 回路能量平衡

由于集中供熱管網(wǎng)為閉合的恒定流系統(tǒng),因此根據(jù)能量守恒定律,沿任一回路方向各管道壓降的代數(shù)和為零。供熱管網(wǎng)中的水力壓降也同樣可以類比電學中的電力壓降,所以沿任意回路方向各個管道壓降的代數(shù)和為零。式(2)為M階列向量ΔH,表示系統(tǒng)中各管道的壓降。因此結合基爾霍夫電壓定律,可以得到通用的回路水力能量方程(3)

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

(2)

(3)

同理,對于管網(wǎng)中任一回路也均應滿足溫降之和為零,因此可以構建回路熱力能量方程(5)。式(4)為M階列向量ΔT,表示各個管道的溫降;bsj表示管道與回路的關聯(lián)信息,為基本回路矩陣B中的元素

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

(4)

(5)

1.2 熱力模型

1.2.1 熱力特性

蒸汽管道可近似看作單層圓筒壁,而實際的蒸汽管道是由n層不同材料組成的,因此多層圓筒壁的總熱阻等于各層熱阻之和。蒸汽管道的總傳熱系數(shù)主要與蒸汽和管道內壁的對流換熱、保溫層的導熱、以及空氣(架空)或土壤(直埋)與管道外護層的換熱等因素有關。由于鋼管管壁導熱熱阻遠小于保溫層導熱熱阻,一般可忽略不計,所以傳熱系數(shù)[11]可按式(6)計算,式中參數(shù)如圖2中定義所示。其中保溫層的導熱主要與保溫材質和保溫層厚度有關系,在此不做詳述,其他可由式(6)-式(10)推導分析。

圖2 保溫管道結構示意圖

(6)

管道內蒸汽在熱源壓力影響下向前受迫流動,對流換熱的實驗關聯(lián)式采用式(7)迪圖斯-貝爾特公式計算。根據(jù)努塞爾數(shù)的定義,蒸汽與管道內壁的對流換熱系數(shù)可由式(8)計算

(7)

(8)

式中Nu——努塞爾數(shù);

Re——雷諾數(shù);

Pr——普朗特數(shù);

hm——蒸汽與管道內壁的對流換熱系數(shù)/W·(m2·K)-1。

環(huán)境與架空管道外護層表面的對流換熱系數(shù),一般可根據(jù)具體風速按式(9)校核。土壤與地埋管道的換熱系數(shù),按福爾赫蓋伊默推導的傳熱學理論公式計算,如式(10)所示

(9)

(10)

式中 h0——管道外護層對流換熱系數(shù)/w·(m2·K)-1;

λt——土壤的導熱系數(shù)/W/(m·K)-1;

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

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

1.2.2 熱力模型與求解

蒸汽各節(jié)點溫度主要取決于管道和附件的沿途溫降,而沿途溫降主要與沿途管道的散熱、蒸汽流量、管道長度以及蒸汽的定壓比熱有關。由此根據(jù)熱平衡原理,結合節(jié)點流量平衡方程、回路能量方程和溫降計算公式,建立熱力求解方程組(11)[12]。如圖3所示,溫度插值策略求解過程中管網(wǎng)熱力計算時的流量和壓力以水力計算的結果為基礎,對節(jié)點溫度反復修正。不僅可以減小復雜水力條件下熱力參數(shù)對水力參數(shù)的影響,提高穩(wěn)定性,而且采用兩個小方程組代替一個大方程組,可以加快求解速度。

圖3 溫度插值策略求解流程圖

(11)

式中A0——由aij組成的(N-1)×M降階關聯(lián)矩陣;

F——由管道溫降系數(shù)組成的M階節(jié)點對角矩陣;

AT——矩陣A的轉置矩陣;

T——節(jié)點溫度向量,T=(t1,t2,t3,…,tN-1)T。

方程組(11)共有2M+N-1個方程,未知各管道的流量、溫降以及各節(jié)點的溫度共2M+N-1個,即方程個數(shù)與未知參數(shù)數(shù)目相等,對該方程組求解可得節(jié)點溫度方程(13)

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

(12)

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

(13)

通過上述熱力計算模型求解分析可知,在節(jié)點溫度及管道溫降求解過程中管道流量對計算結果有著較為重要的影響,為使系統(tǒng)計算得到的結果更加符合實際并減小誤差,應首先對系統(tǒng)流量進行分配。

1.3 水力模型

1.3.1 水力特性

供熱管道的壓力損失等于始、末兩斷面間沿程和局部阻力損失之和,如式(14)所示。在單位長度范圍內,假設管道摩擦阻力系數(shù)和管內蒸汽密度不隨管道長度變化,沿程阻力損失可由Darcy公式(15)計算

ΔHh=ΔHl+ΔHf

(14)

(15)

式中 ΔHh——壓力損失/Pa;

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

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

L——長度/m;

λ——阻力系數(shù);

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

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

大部分蒸汽的流動都屬于紊流過渡區(qū),如式(16)所示,Colebrook-White公式給出了管道當量粗糙度、雷諾數(shù)、阻力系數(shù)的關系,適用于工業(yè)管道的三個阻力區(qū),又稱為紊流λ的綜合公式。阻力系數(shù)代入Darcy公式,可以得到目前理論上最精確的壓力損失計算結果。蒸汽流經彎頭、三通等部件所產生的局部阻力損失,通常采用局部阻力系數(shù)和當量長度近似方法計算。結合管道內蒸汽流速計算公式(17),可以得到局部阻力損失計算公式(18)

(16)

(17)

(18)

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

Le——當量長度/m。

假定管道內的凝結水均在節(jié)點處生成排出,蒸汽流量不變,管道內流體的壓降和流量的關系服從二次冪規(guī)律。蒸汽管網(wǎng)穩(wěn)態(tài)水力計算壓力損失數(shù)學表達式可以歸納為式(19)

(19)

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

1.3.2 水力模型與求解

對于已知的管網(wǎng)系統(tǒng),參照電路理論中的基爾霍夫定律,結合阻力損失公式、節(jié)點流量平衡和回路能量平衡方程等建立通用的蒸汽管網(wǎng)流量分配方程組(20)[13]

(20)

式中S=diag(s1,s2,s3,…,sM)——管道阻力特性系數(shù)對角矩陣;

H=(H1,H2,H3,…HN-1)T——節(jié)點壓力向量;

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

由方程組(20)可知該管道流量分配模型中管道流量滿足節(jié)點流量平衡方程,同時由方程組(11)可知熱力計算過程中管道流量同樣滿足節(jié)點流量平衡方程。因此,熱力計算時的管道流量采用該模型計算得出的結果的方法可行。式(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階管網(wǎng)系數(shù)矩陣。根據(jù)節(jié)點分類的不同,劃分管網(wǎng)系數(shù)矩陣,可以得到控制性方程組(23),目前被認為是最簡潔、高效的管網(wǎng)計算模型[14]

(23)

式中A11、0、A21、A22——系數(shù)矩陣劃分后的分塊矩陣;

A10=[M;N-n0]——管道與已知壓力節(jié)點的連接矩陣;

A21——阻力特性系數(shù)對角矩陣S;

QT=[q1,q2,…,qN-n0]——已知節(jié)點需求量;

HT=[H1,H2,…,Hn0]——未知節(jié)點壓力。

為方便求解,將非線性方程進行線性化處理,用求解線性方程組的方法逐步逼近非線性方程組的解,并利用迭代的方法消除誤差。具體數(shù)值求解過程為式(24)～式(26)。利用Newton-Jacobi迭代方法可求取關于管道流量與節(jié)點壓力的修正值,基于此方程,系統(tǒng)模型可由式(26)逐步迭代求解

(24)

(25)

(26)

以上為管網(wǎng)系統(tǒng)水力計算的主要過程,將管道壓力降與管道流量之間的非線性關系轉化為線性關系,迭代修正參數(shù)與基值之間的誤差,使前后兩次的值不斷接近,獲得系統(tǒng)的壓力和流量分布。針對熱力計算,可以采用相似的數(shù)值求解方法。

2.2 案例分析

為驗證本文所建立模型的有效性,對某工業(yè)園區(qū)進行仿真分析。根據(jù)供熱系統(tǒng)結構簡化拓撲關系,如圖4所示:由29個節(jié)點、28條管道組成,每條管道末端為熱用戶,主要有16家用熱企業(yè)。該蒸汽管網(wǎng)為枝狀,全長約27 km,管徑范圍從DN70到 DN600。熱源廠供給壓力為1.6 MPa,溫度為250 ℃。管網(wǎng)具體設計參數(shù)如表1所示。

圖4 園區(qū)供熱管網(wǎng)規(guī)劃示意圖

圖5 節(jié)點運行和仿真數(shù)據(jù)對比

圖6 壓力和溫度誤差分析

表1 園區(qū)管網(wǎng)設計參數(shù)

由于該工業(yè)園區(qū)尚未對凝結水量進行監(jiān)測,提供流量參數(shù)僅有熱源供給量和熱用戶需求量,所建立的一維穩(wěn)態(tài)模型假設蒸汽過熱輸送且為單一均相,無法對管道凝結水量進行評估。因此不再對管道蒸汽和凝結水流量討論,仿真結果如表2所示。

表2 管網(wǎng)系統(tǒng)仿真結果

續(xù)表2

通過以上圖表分析可知,在大部分節(jié)點處,仿真計算結果與實際運行數(shù)據(jù)吻合度較好,溫度和壓力誤差大多在5%以內。節(jié)點壓力和溫度,在上游靠近熱源點處與實際運行數(shù)據(jù)較為接近,但隨著與熱源點距離的增加,相對誤差不斷增大,節(jié)點F2、G、f2、F、f的誤差相對較大。其中F點壓力誤差最大為 5.74%,a點壓力誤差最小為 0.193%,平均壓力誤差為 3.05%;f點溫度誤差最大為 5.45%,C11點溫度誤差最小為 1.10%,平均溫度誤差為2.69%。這是由于節(jié)點參數(shù)是通過與熱源點的差值迭代求得,因此隨著距離的增加計算誤差不斷積累變大。

3 結論

本文通過對蒸汽管網(wǎng)壓力損失特性的分析,結合蒸汽管道與外界環(huán)境的傳熱過程特征和熱平衡原理,得到了描述蒸汽管網(wǎng)水力和熱力計算的數(shù)學模型。實現(xiàn)了蒸汽傳輸過程的完整數(shù)值描述和整體線性化求解。

(1)線性化求解過程中,首先初始化系統(tǒng)輸入狀態(tài)參數(shù),通過水力模型求解壓力和流量,其次將計算結果傳遞給熱力模型求解溫度,然后將參數(shù)值反饋給水力模型再次求解,如此循環(huán)至仿真精度為止。在對數(shù)值精度影響較小的前提下,保證管網(wǎng)的快速求解和收斂性。

(2)通過節(jié)點質量平衡和回路能量平衡關系,將模型轉化為以節(jié)點壓力或溫度為變量的方程進行迭代求解,得到了管網(wǎng)系統(tǒng)運行參數(shù)。通過實例驗證,大部分仿真誤差在5%以內,基本滿足工程應用要求?？紤]模型的簡化處理、蒸汽的可壓縮性、管件老化等對計算結果產生的影響,該模型可為管網(wǎng)系統(tǒng)的優(yōu)化運行提供理論依據(jù)。