王海超, 謝子蟬, Risto Lahdelma, 華鵬敏, 周 揚(yáng)

(1.大連理工大學(xué)土木工程學(xué)院,遼寧大連116024;2.阿爾托大學(xué)科學(xué)學(xué)院, 芬蘭艾斯堡11100)

1 概述

供熱系統(tǒng)的動(dòng)態(tài)仿真問題受到了很多的關(guān)注,已開展研究也很多。因?yàn)槠洳粌H為供熱系統(tǒng)的規(guī)劃設(shè)計(jì)及運(yùn)行調(diào)度服務(wù),也可以為當(dāng)前的智慧供熱系統(tǒng)及人工智能算法引入供熱智能調(diào)控服務(wù),為實(shí)現(xiàn)更加精準(zhǔn)的供熱調(diào)度提供了基礎(chǔ)。

動(dòng)態(tài)熱力仿真中,溫度反饋是實(shí)現(xiàn)閉環(huán)控制的關(guān)鍵[1],動(dòng)態(tài)熱力仿真是獲得溫度反饋的一種經(jīng)濟(jì)可行的方法[2-3]。熱網(wǎng)模型應(yīng)該準(zhǔn)確模擬熱傳播,這對(duì)負(fù)荷調(diào)節(jié)及運(yùn)行很重要[4-6]。熱在管網(wǎng)中的傳播受傳輸延時(shí)和熱損失兩個(gè)主要參數(shù)的影響。從機(jī)理上看,熱傳播的速度略滯后于流動(dòng)速度,但不會(huì)相差太多[7]。在典型的運(yùn)行工況下,傳輸時(shí)間延遲為min級(jí)和h級(jí)[8]。根據(jù)熱網(wǎng)結(jié)構(gòu)、規(guī)模和運(yùn)行條件不同,傳輸中的熱損失從年供熱量的5%到20%不等[2]。因此,應(yīng)該準(zhǔn)確地計(jì)算每條管道的熱傳播。

在流體力學(xué)中,描述質(zhì)點(diǎn)運(yùn)動(dòng)的方法有兩種:歐拉法和拉格朗日法。在歐拉法中,方程建立在固定的位置上,網(wǎng)格與質(zhì)點(diǎn)的速度無關(guān)。與此相比,拉格朗日法跟蹤運(yùn)動(dòng)的控制體粒子,對(duì)于動(dòng)態(tài)熱力仿真可能更有效。歐拉法由于具有求解偏微分方程組的固定控制量的優(yōu)點(diǎn),在數(shù)值研究中得到了廣泛應(yīng)用。Wang等人[6]提出了一階隱式迎風(fēng)模型,將其與基于特征線法的模型[9]進(jìn)行了比較。結(jié)果表明,特征線法的計(jì)算速度是一階隱式迎風(fēng)模型的兩倍,但精度較差。為了進(jìn)一步提高精度,Zheng等人[10]提出了一種新的測(cè)量方法,介紹了二階隱式迎風(fēng)方法和半隱式快速方法,并與一階方法進(jìn)行了比較。3個(gè)模型都使用了相同的實(shí)驗(yàn)數(shù)據(jù),這表明二階隱式迎風(fēng)方法是最有效的方法,出口溫度誤差在±1.0 ℃以內(nèi)。Betancourt等人[3]提出了一種改進(jìn)的有限體積法(Finite Volume Method,FVM),該方法計(jì)算了一個(gè)有6個(gè)節(jié)點(diǎn)的網(wǎng)絡(luò)中的節(jié)點(diǎn)溫度。結(jié)果表明,該方法的計(jì)算速度是傳統(tǒng)有限體積法的3.5倍,但精度與Courant準(zhǔn)則有關(guān)。Courant準(zhǔn)則是數(shù)值收斂的必要條件,它要求時(shí)間步長小于波傳播到相鄰網(wǎng)格的持續(xù)時(shí)間[11]。

歐拉法計(jì)算時(shí)長與迭代步數(shù)和每一步的復(fù)雜度有關(guān),精度與時(shí)間和空間的離散化密切相關(guān)。當(dāng)采用粗網(wǎng)格時(shí),精度表現(xiàn)不佳,在不犧牲精度的情況下,需要大幅提高求解偏微分方程組的速度。還有研究者深入研究了時(shí)間步長和空間步長尺度的適當(dāng)設(shè)置,以平衡精度和模擬時(shí)間[6,10,12],但由于沒有通用的模型設(shè)置規(guī)則,該方法的大規(guī)模應(yīng)用有一定的阻礙。

拉格朗日法的一個(gè)著名的簡(jiǎn)化模型稱為節(jié)點(diǎn)法,由Benonysson于1991年首次提出[13]。其想法是跟蹤控制體積從管道的起點(diǎn)到終點(diǎn)的行程時(shí)間。當(dāng)水進(jìn)入管道時(shí),會(huì)產(chǎn)生可變控制量。這意味著空間離散化由每個(gè)時(shí)間步長的水的體積流量決定。進(jìn)水口平均加權(quán)水溫在控制體積離開管道之前保持不變。該時(shí)間步長中的節(jié)點(diǎn)溫度由更新后的控制體的溫度給出。后來發(fā)展了節(jié)點(diǎn)法和解析解的組合模型,該模型在精度和計(jì)算速度上具有相似的性能[14]。比如TERMIS軟件是基于節(jié)點(diǎn)法開發(fā)的,但沒有考慮管材的熱容[8]。在Gabrielaitiene等人的研究[15-16]中,對(duì)節(jié)點(diǎn)法和TERMIS軟件在實(shí)際項(xiàng)目中的性能進(jìn)行了詳細(xì)的比較。類似的工作可以在基于Modelica的框架中找到[17]。由于空間離散的靈活性,這些基于拉格朗日有限體積的模型比基于歐拉方法的模型更加靈活和準(zhǔn)確[18]。

上述研究表明,在以前的研究中,往往采用統(tǒng)一的時(shí)間步長,但是仍然需要更小的時(shí)間步長來精確地模擬急劇的溫度變化,這會(huì)導(dǎo)致計(jì)算速度降低,還會(huì)導(dǎo)致數(shù)值不穩(wěn)定[17]。當(dāng)前如何應(yīng)用可變的時(shí)間步長和空間步長來實(shí)現(xiàn)快速、準(zhǔn)確的動(dòng)態(tài)熱力仿真是一個(gè)急需解決的問題。為此,本文提出了基于離散事件模擬(Discrete Event Simulation,DES)[19]的供熱管道快速動(dòng)態(tài)熱力仿真方法,這是首次將DES應(yīng)用于供熱管道動(dòng)態(tài)仿真的研究。該方法考慮了水在管道中的熱損失和時(shí)滯。它可以高效、準(zhǔn)確地處理不同的入口溫度和流量,包括零流量。本文對(duì)管道內(nèi)的溫度傳播進(jìn)行了詳細(xì)的計(jì)算,利用測(cè)量數(shù)據(jù)對(duì)模型進(jìn)行了驗(yàn)證。

2 供熱管道離散事件仿真模型建立

溫度傳播計(jì)算基于活塞流模型,假定管道的所有橫截面速度均勻。為了快速計(jì)算,忽略了溫度邊界和水力邊界的差異。對(duì)于熱傳遞的影響,假定環(huán)境溫度不變計(jì)算熱損失。不考慮管壁與水之間的摩擦產(chǎn)生熱量,不考慮軸向傳熱。

2.1 離散事件仿真思想

離散事件仿真的核心思想是將供熱管網(wǎng)中的并行過程模擬為一系列事件。水在管內(nèi)的一個(gè)流動(dòng)狀態(tài)變化視為一個(gè)事件,因此供熱管網(wǎng)運(yùn)行是多個(gè)連續(xù)事件的組合。事件可以被創(chuàng)建、修改,甚至刪除,即事件序列(Event Quene,EQ)可根據(jù)事件驅(qū)動(dòng)閾值更新,與固定時(shí)間步長的仿真方法不同,DES中的仿真時(shí)間是從當(dāng)前事件時(shí)間跳躍到下一個(gè)事件時(shí)間,這意味著時(shí)間步長取決于事件的激活,因此是可變長度的。

該模擬基于跟蹤水力前沿在管道中的運(yùn)動(dòng)。水邊界是在假定為柱塞流的情況下在管道中流動(dòng)的無限短的水段。DES模型包括3種類型的事件:入口溫度變化事件、流量變化事件(也稱流速變化事件)、水邊界到達(dá)管道末端事件(簡(jiǎn)稱水邊界到達(dá)事件)。入口溫度變化事件和流量變化事件可以獨(dú)立地發(fā)生,由輸入數(shù)據(jù)給出,在模擬過程中將生成水邊界到達(dá)管道末端事件?？梢娫趦蓚€(gè)流量變化事件之間,流量保持不變。

當(dāng)入口溫度變化或流速變化時(shí),在管道的入口處形成水邊界。水邊界穿過管道,直到到達(dá)管道末端,這構(gòu)成了到達(dá)事件。水邊界的順序由先進(jìn)先出(First in First out,FIFO)隊(duì)列管理。最先到達(dá)管道末端的水邊界稱為第一邊界,而最新創(chuàng)造的水邊界是最后的邊界。開始的時(shí)候隊(duì)列是空的,并且管道被初始化為基于恒定流量和入口溫度的穩(wěn)態(tài)操作。當(dāng)FIFO隊(duì)列非空時(shí),第一邊界的到達(dá)事件始終被調(diào)度在事件隊(duì)列中。當(dāng)邊界到達(dá)第一位置時(shí),計(jì)算到達(dá)事件的激活時(shí)間,并在流速變化時(shí)更新該時(shí)間。

DES仿真的核心思想見圖1。圖1中t表示時(shí)間,共包含4個(gè)事件:t1處1個(gè)入口溫度變化事件,t2和t3處的2個(gè)流速變化事件,t4處的1個(gè)水邊界到達(dá)事件。圖1中,F1、F2、F3分別表示不同時(shí)間創(chuàng)建的水力前沿。所謂水力前沿,指熱水在管道中溫度或流速發(fā)生變化的位置,隨著時(shí)間推移,水力前沿會(huì)向管道末端流動(dòng)。圓柱體表示管道。

圖1 DES仿真的核心思想(軟件截圖)

在t1時(shí)刻,入口溫度改變,創(chuàng)建第1個(gè)水力前沿F1,并根據(jù)當(dāng)前流速計(jì)算相應(yīng)的到達(dá)事件的時(shí)間。在t2時(shí)刻,流速改變,創(chuàng)建第2個(gè)水力前沿F2,并根據(jù)新的流速重新計(jì)算F1的到達(dá)事件時(shí)間。在t3時(shí)刻,再次改變流速,創(chuàng)建水力前沿F3,并重新計(jì)算F1的到達(dá)事件時(shí)間。在t4時(shí)刻,F1到達(dá)管道末端,并從FIFO隊(duì)列中移除,F2變成第1個(gè)前沿,再根據(jù)當(dāng)前流速,計(jì)算其到達(dá)事件時(shí)間,以此類推。

2.2 行程時(shí)間計(jì)算

水邊界的傳輸時(shí)間取決于它們?cè)谕ㄟ^管道時(shí)所經(jīng)歷的可變流速集。根據(jù)管道的橫截面面積和水密度,根據(jù)質(zhì)量流量計(jì)算給定時(shí)刻的流速,計(jì)算式為:

(1)

式中v——流速,m/s

qm——水的質(zhì)量流量,kg/s

ρ——水密度,kg/m3

A——管道橫截面面積,m2

每個(gè)水邊界的行程時(shí)間是水邊界到達(dá)管道末端時(shí)間(簡(jiǎn)稱到達(dá)時(shí)間)和水邊界創(chuàng)造時(shí)間之差值,即:

ttra=tout-tin

(2)

式中ttra——水邊界行程時(shí)間,s

tout——水邊界到達(dá)時(shí)間,s

tin——水邊界創(chuàng)造時(shí)間,s

對(duì)1個(gè)水邊界而言,第1個(gè)水邊界的到達(dá)總是按該時(shí)刻計(jì)算的預(yù)期到達(dá)時(shí)間被安排在事件隊(duì)列中。假設(shè)流速v>0,則第i個(gè)水邊界的預(yù)期到達(dá)時(shí)間計(jì)算如下:

(3)

式中tiout——第i個(gè)水邊界的到達(dá)時(shí)間,s

tiin——第i個(gè)水邊界的創(chuàng)造時(shí)間,s

di——第i個(gè)水邊界到前(i-1)個(gè)水邊界的剩余行程距離,m

如果不存在先前的水邊界,則將di設(shè)置為管道長度L。在前1個(gè)水邊界到達(dá)之后,或者在創(chuàng)建水邊界時(shí)FIFO隊(duì)列為空的情況下,下1個(gè)水邊界成為第1個(gè)。當(dāng)?shù)趇個(gè)水邊界成為第1個(gè)時(shí),它到管道末端(出口)的剩余距離等于di。

當(dāng)流速變化時(shí),需要更新水邊界到達(dá)時(shí)間。新的水邊界到達(dá)時(shí)間tnewout計(jì)算式為:

(4)

式中tnewout——新的水邊界到達(dá)時(shí)間,s

tc——流速變化時(shí)刻,s

vnew——新的流速,m/s

在時(shí)間t到達(dá)出口的控制體粒子行程時(shí)間ttra(t)是一個(gè)分段線性函數(shù),且只要流量為正即流向不變化,該函數(shù)就是連續(xù)的。

這意味著控制體粒子從tin到新的創(chuàng)造時(shí)間tnewin的歷史入口速度的積分與控制體粒子從tout到tnewout的到達(dá)速度的積分相同。由于兩個(gè)相鄰水邊界之間的控制體粒子入口速度是恒定的,如第2.1節(jié)所述,當(dāng)流速在tout和tnewout之間保持不變時(shí),得到計(jì)算式(5):

(5)

式中tnewin——新的創(chuàng)造時(shí)間,s

vout——控制體粒子離開流速,m/s

vin——控制體粒子進(jìn)入流速,m/s

新的行程時(shí)間tnewtra可以根據(jù)水邊界行程時(shí)間ttra計(jì)算,計(jì)算式見式(6)。雖然式(6)是針對(duì)正流速推導(dǎo)的,但它也適用于零流速(vout= 0)。式(6)如下:

(6)

式中tnewtra——新的行程時(shí)間,s

Δt——新的水邊界到達(dá)時(shí)間與原到達(dá)時(shí)間的差,s

在水邊界到達(dá)時(shí)間以及在流速變化時(shí),更新行程時(shí)間和斜率。因此,當(dāng)兩個(gè)到達(dá)事件之間的流速保持恒定時(shí),兩個(gè)到達(dá)事件之間的行程時(shí)間是線性函數(shù)。流速的變化使行程時(shí)間在整個(gè)模擬時(shí)間內(nèi)成為分段線性函數(shù)。當(dāng)流速為正值時(shí),行程時(shí)間函數(shù)是連續(xù)的。零流速會(huì)在行程時(shí)間函數(shù)中產(chǎn)生不連續(xù)的跳躍。

2.3 溫度計(jì)算

在水密度和比定壓熱容不變的情況下,能量平衡方程見式(7):

(7)

式中cp——水的比定壓熱容,J/(kg·K)

t——時(shí)間,s

x——位置,m

θ(t,x)——位置x處控制體粒子在時(shí)間t時(shí)的溫度,℃

q(t,x)——位置x處控制體粒子在時(shí)間t時(shí)的從水到環(huán)境的熱損失率,W

式(7)左邊是內(nèi)部熱能關(guān)于時(shí)間t和位置x的導(dǎo)數(shù)項(xiàng),右邊是熱損失率q(t,x)。q(t,x)為正表示從水到環(huán)境熱傳遞。

本研究在熱損失計(jì)算中采用了水到環(huán)境的管道總熱阻R,避免了對(duì)不同層數(shù)的迭代計(jì)算。式(7)可簡(jiǎn)化為:

(8)

或:

(9)

式中θ——控制體粒子溫度,℃

θamb——環(huán)境溫度,可以是環(huán)境空氣溫度,也可以是地面溫度,℃

R——管道總熱阻,m·K/W

對(duì)式(9)積分可得到出口溫度方程,為:

(10)

式中θout——控制體粒子在行程時(shí)間后到達(dá)出口時(shí)的當(dāng)前溫度,℃

θin——控制體粒子初始溫度,℃

因此,管道中的控制體粒子溫度是漸近接近環(huán)境溫度的時(shí)間的指數(shù)函數(shù),單個(gè)控制體粒子通過管道時(shí)的溫度分布見圖2。為了便于比較,采用線性插值法計(jì)算兩個(gè)相鄰斷點(diǎn)之間的溫度。當(dāng)指數(shù)函數(shù)接近線性時(shí),即當(dāng)溫降較小時(shí),線性近似是合理的。

圖2 單個(gè)控制體粒子通過管道時(shí)的溫度分布

采用漸變模式來表示時(shí)間區(qū)間內(nèi)的入口溫度分布。即入口溫度在每個(gè)時(shí)間區(qū)間內(nèi)逐漸變化。采用指數(shù)插值法時(shí),入口溫度分布為分段指數(shù)函數(shù),而采用線性插值法時(shí),入口溫度分布為分段線性函數(shù)。

2.4 熱損失計(jì)算

熱損失可直接由控制體粒子溫降計(jì)算,也可間接由能量平衡計(jì)算。本文選擇能量平衡法,其中參考溫度為環(huán)境溫度。在該模型中,可以計(jì)算任意時(shí)間段的熱量損失。

在開始計(jì)算時(shí)間t0和結(jié)束計(jì)算時(shí)間tend之間,由環(huán)境造成的水的熱損失見式(11)～(14):

(11)

(12)

(13)

(14)

式中Qloss——在t0和tend之間管長為L的管道的熱損失,J

t0——開始計(jì)算時(shí)間,s

tend——結(jié)束計(jì)算時(shí)間,s

L——管道長度,m

Qin——入口能量,J

Qout——出口能量,J

U(tend)——tend時(shí)刻熱力學(xué)能,J

U(t0)——t0時(shí)刻熱力學(xué)能,J

v(t) ——t時(shí)刻流速,m/s

θin(t)——t時(shí)刻控制體粒子初始溫度,℃

θout(t)——t時(shí)刻控制體粒子在行程時(shí)間后到達(dá)出口時(shí)的當(dāng)前溫度,℃

U(t)——t時(shí)刻熱力學(xué)能,J

式(11)的右側(cè)項(xiàng)表示入口能量、出口能量和內(nèi)部能量的變化。入口能量Qin和出口能量Qout是累積變量,在模擬過程中會(huì)進(jìn)行更新。Qin和Qout可以分別更新。在產(chǎn)生新的水邊界時(shí)更新Qin,而在出口溫度函數(shù)的斷點(diǎn)處更新Qout。t時(shí)刻熱力學(xué)能U(t)是一個(gè)狀態(tài)變量,它只在t0和tend處計(jì)算。

由于斷點(diǎn)之間是恒定流速,式(12)、(13)中的積分是分段指數(shù)函數(shù),其斷點(diǎn)與溫度函數(shù)的斷點(diǎn)匹配。因此,我們將Qin和Qout的集成分解為幾個(gè)片段。同樣,在任意時(shí)刻t,管道內(nèi)的溫度分布是距管道入口的距離x的分段指數(shù)函數(shù)。

在線性插值法中,假設(shè)入口溫度是t的分段線性函數(shù),出口溫度分布用分段線性函數(shù)逼近。因?yàn)閮蓚€(gè)斷點(diǎn)之間的流量恒定,所以積分在式(12)、(13)中由于線性被積函數(shù)的緣故可以被簡(jiǎn)化。為了計(jì)算熱力學(xué)能,需要從入口到出口的水溫分布。作為位置x的函數(shù)的水溫分布也由分段線性函數(shù)來近似。

上述DES模型稱為漸變線性近似模型,即GL模型,也是可以建立的DES模型中的基本模型。

3 模型驗(yàn)證

本研究利用石家莊一條供熱管道的24 h入口溫度、出口溫度和流量數(shù)據(jù)對(duì)GL模型進(jìn)行驗(yàn)證[14]。

供熱管道入口溫度在88.4～97.9 ℃范圍內(nèi)變化,流量在9 012.4～9 761.8 m3/h范圍內(nèi)變化,測(cè)量時(shí)間步長為5 min,溫度傳感器的分度值為0.4 ℃,流量計(jì)分度值為0.01 m3/h[19]。根據(jù)提取的數(shù)據(jù),有256個(gè)入口溫度變化事件和258個(gè)流量變化事件作為模型的輸入。

案例供熱管道入口流量數(shù)據(jù)[15]見圖3,案例供熱管道其他參數(shù)[15]見表1。管道的初始狀態(tài)未知,但假設(shè)首次測(cè)量的入口溫度為88.5 ℃、流量為9 761.8 m3/h。

圖3 案例供熱管道流量數(shù)據(jù)(軟件截圖)

表1 案例供熱管道其他參數(shù)

① 出口溫度誤差

模擬過程中前1.58 h為初始化階段,直至第一個(gè)水邊界到達(dá)管道末端。出口溫度測(cè)量值與GL模型模擬值的誤差見圖4。

圖4 出口溫度測(cè)量值與GL模型模擬值的誤差(軟件截圖)

由圖4可以看出,除初始化階段外,出口溫度模擬值和測(cè)量值有較好的一致性。GL模型的最大誤差為±0.6 ℃?？傮w上看,出口溫度模擬值略低于測(cè)量值,這可能是對(duì)熱損失的輕微高估造成的,后續(xù)將進(jìn)一步完善熱損失模型。這與特征線法和隱式迎風(fēng)法[15]相當(dāng),優(yōu)于一階隱式迎風(fēng)法[19],但計(jì)算速度更快。

② 管道水溫分布

GL模型模擬可以提供管道沿線的溫度分布,獲得任意時(shí)刻沿管道的溫度分布曲線。GL模型模擬的9:00和9:30的供熱管道溫度分布見圖5。

圖5 GL模型模擬的9:00和 9:30的供熱管道溫度分布(軟件截圖)

由圖5可以看出,在30 min內(nèi)水力前沿前進(jìn)了約3 053 m,并略有降溫,但溫度降幅很小,不到0.1 °C。

當(dāng)水邊界對(duì)管道溫度分布影響很小時(shí)刪除水邊界是可能的,這涉及如何確定一個(gè)閾值,從而確定是否產(chǎn)生了新的事件。在大型模型中,刪除這種多余的影響不大的水邊界可能會(huì)顯著加快計(jì)算速度,刪除多余水邊界的標(biāo)準(zhǔn)是未來需要進(jìn)一步研究的問題,這涉及計(jì)算量和精度的平衡。

4 結(jié)論及展望

① 提出了一種供熱管道動(dòng)態(tài)熱力仿真方法。該方法基于拉格朗日方法,采用離散事件模擬(Discrete Event Simulation,DES)跟蹤沿管道傳輸?shù)乃η把?是首次將DES應(yīng)用于供熱管道動(dòng)態(tài)仿真的研究。該方法將水在管內(nèi)的1個(gè)流動(dòng)狀態(tài)變化視為1個(gè)事件,將整個(gè)供熱管網(wǎng)模擬分為多個(gè)連續(xù)事件的組合。事件可以被創(chuàng)建、修改、刪除。仿真的時(shí)間步長和空間步長可變,在達(dá)到相同計(jì)算精度的前提下,計(jì)算速度可大幅提升。

② 對(duì)GL模型進(jìn)行案例驗(yàn)證,結(jié)果表明該方法完全可行。除初始化階段外,出口溫度模擬值和測(cè)量值有較好的一致性。GL模型出口溫度模擬值的最大誤差為±0.6 ℃,但出口溫度模擬值略低于測(cè)量值,這可能是對(duì)熱損失的輕微高估造成的,后續(xù)將進(jìn)一步完善熱損失模型。

③ GL模型可以隨時(shí)提供管道沿線的溫度分布,獲得任意時(shí)刻沿管道的溫度分布曲線。

④ 通過研究同時(shí)發(fā)現(xiàn),當(dāng)水邊界對(duì)管道溫度分布影響很小時(shí)刪除水邊界是可能的,這涉及如何確定一個(gè)閾值,從而確定是否產(chǎn)生了新的事件。在大型模型中,刪除這種多余的影響不大的水邊界可能會(huì)顯著加快計(jì)算速度,刪除多余水邊界的標(biāo)準(zhǔn)是未來需要進(jìn)一步研究的問題,這涉及計(jì)算量和精度的平衡。因此后續(xù)將繼續(xù)開發(fā)基于該核心思想的平臺(tái)產(chǎn)品,并研究該方法在大規(guī)模熱網(wǎng)和環(huán)狀熱網(wǎng)中的適用性和仿真速度。