基于溫度插值求解的管網仿真研究

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

(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%。這是由于節點參數是通過與熱源點的差值迭代求得,因此隨著距離的增加計算誤差不斷積累變大。

3 結論

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

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

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