基于溫度插值求解的管網(wǎng)仿真研究

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

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

0 引言

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

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

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

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

1 數(shù)學(xué)模型

1.1 質(zhì)能平衡

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

1.1.1 節(jié)點(diǎn)質(zhì)量平衡

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

圖1 節(jié)點(diǎn)質(zhì)量平衡圖

(1)

式中A——流出或流入對(duì)應(yīng)關(guān)聯(lián)矩陣;

aij——A中的元素;

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

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

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

1.1.2 回路能量平衡

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

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

(2)

(3)

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

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

(4)

(5)

1.2 熱力模型

1.2.1 熱力特性

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

圖2 保溫管道結(jié)構(gòu)示意圖

(6)

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

(7)

(8)

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

Re——雷諾數(shù);

Pr——普朗特?cái)?shù);

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

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

(9)

(10)

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

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

w——風(fēng)速/m·s-1;

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

1.2.2 熱力模型與求解

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

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

(11)

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

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

AT——矩陣A的轉(zhuǎn)置矩陣;

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

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

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

(12)

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

(13)

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

1.3 水力模型

1.3.1 水力特性

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

ΔHh=ΔHl+ΔHf

(14)

(15)

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

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

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

L——長(zhǎng)度/m;

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

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

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

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

(16)

(17)

(18)

式中ε——管道相對(duì)粗糙度/m;

Le——當(dāng)量長(zhǎng)度/m。

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

(19)

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

1.3.2 水力模型與求解

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

(20)

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

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

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

由方程組(20)可知該管道流量分配模型中管道流量滿足節(jié)點(diǎn)流量平衡方程,同時(shí)由方程組(11)可知熱力計(jì)算過(guò)程中管道流量同樣滿足節(jié)點(diǎn)流量平衡方程。因此,熱力計(jì)算時(shí)的管道流量采用該模型計(jì)算得出的結(jié)果的方法可行。式(21)和(22)為方程組(21)化簡(jiǎn)求解過(guò)程

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

(21)

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

(22)

2 模型線性化求解與案例分析

2.1 模型線性化求解

上述模型僅包含流量、壓力和溫度三個(gè)基礎(chǔ)變量,式(12)和(21)中(A·(F|G|)-1·AT)-1、(A·(S|G|)-1·AT)-1實(shí)際為N·N階管網(wǎng)系數(shù)矩陣。根據(jù)節(jié)點(diǎn)分類(lèi)的不同,劃分管網(wǎng)系數(shù)矩陣,可以得到控制性方程組(23),目前被認(rèn)為是最簡(jiǎn)潔、高效的管網(wǎng)計(jì)算模型[14]

(23)

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

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

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

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

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

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

(24)

(25)

(26)

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

2.2 案例分析

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

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

圖5 節(jié)點(diǎn)運(yùn)行和仿真數(shù)據(jù)對(duì)比

圖6 壓力和溫度誤差分析

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

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

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

續(xù)表2

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

3 結(jié)論

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

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

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