基于離散事件模擬的熱網快速動態熱力仿真

王海超, 謝子蟬, Risto Lahdelma, 華鵬敏, 周 揚

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

1 概述

供熱系統的動態仿真問題受到了很多的關注,已開展研究也很多。因為其不僅為供熱系統的規劃設計及運行調度服務,也可以為當前的智慧供熱系統及人工智能算法引入供熱智能調控服務,為實現更加精準的供熱調度提供了基礎。

動態熱力仿真中,溫度反饋是實現閉環控制的關鍵[1],動態熱力仿真是獲得溫度反饋的一種經濟可行的方法[2-3]。熱網模型應該準確模擬熱傳播,這對負荷調節及運行很重要[4-6]。熱在管網中的傳播受傳輸延時和熱損失兩個主要參數的影響。從機理上看,熱傳播的速度略滯后于流動速度,但不會相差太多[7]。在典型的運行工況下,傳輸時間延遲為min級和h級[8]。根據熱網結構、規模和運行條件不同,傳輸中的熱損失從年供熱量的5%到20%不等[2]。因此,應該準確地計算每條管道的熱傳播。

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

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

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

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

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

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

2.1 離散事件仿真思想

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

該模擬基于跟蹤水力前沿在管道中的運動。水邊界是在假定為柱塞流的情況下在管道中流動的無限短的水段。DES模型包括3種類型的事件:入口溫度變化事件、流量變化事件(也稱流速變化事件)、水邊界到達管道末端事件(簡稱水邊界到達事件)。入口溫度變化事件和流量變化事件可以獨立地發生,由輸入數據給出,在模擬過程中將生成水邊界到達管道末端事件。可見在兩個流量變化事件之間,流量保持不變。

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

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

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

在t1時刻,入口溫度改變,創建第1個水力前沿F1,并根據當前流速計算相應的到達事件的時間。在t2時刻,流速改變,創建第2個水力前沿F2,并根據新的流速重新計算F1的到達事件時間。在t3時刻,再次改變流速,創建水力前沿F3,并重新計算F1的到達事件時間。在t4時刻,F1到達管道末端,并從FIFO隊列中移除,F2變成第1個前沿,再根據當前流速,計算其到達事件時間,以此類推。

2.2 行程時間計算

水邊界的傳輸時間取決于它們在通過管道時所經歷的可變流速集。根據管道的橫截面面積和水密度,根據質量流量計算給定時刻的流速,計算式為:

(1)

式中v——流速,m/s

qm——水的質量流量,kg/s

ρ——水密度,kg/m3

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

每個水邊界的行程時間是水邊界到達管道末端時間(簡稱到達時間)和水邊界創造時間之差值,即:

ttra=tout-tin

(2)

式中ttra——水邊界行程時間,s

tout——水邊界到達時間,s

tin——水邊界創造時間,s

對1個水邊界而言,第1個水邊界的到達總是按該時刻計算的預期到達時間被安排在事件隊列中。假設流速v>0,則第i個水邊界的預期到達時間計算如下:

(3)

式中tiout——第i個水邊界的到達時間,s

tiin——第i個水邊界的創造時間,s

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

如果不存在先前的水邊界,則將di設置為管道長度L。在前1個水邊界到達之后,或者在創建水邊界時FIFO隊列為空的情況下,下1個水邊界成為第1個。當第i個水邊界成為第1個時,它到管道末端(出口)的剩余距離等于di。

當流速變化時,需要更新水邊界到達時間。新的水邊界到達時間tnewout計算式為:

(4)

式中tnewout——新的水邊界到達時間,s

tc——流速變化時刻,s

vnew——新的流速,m/s

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

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

(5)

式中tnewin——新的創造時間,s

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

vin——控制體粒子進入流速,m/s

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

(6)

式中tnewtra——新的行程時間,s

Δt——新的水邊界到達時間與原到達時間的差,s

在水邊界到達時間以及在流速變化時,更新行程時間和斜率。因此,當兩個到達事件之間的流速保持恒定時,兩個到達事件之間的行程時間是線性函數。流速的變化使行程時間在整個模擬時間內成為分段線性函數。當流速為正值時,行程時間函數是連續的。零流速會在行程時間函數中產生不連續的跳躍。

2.3 溫度計算

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

(7)

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

t——時間,s

x——位置,m

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

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

式(7)左邊是內部熱能關于時間t和位置x的導數項,右邊是熱損失率q(t,x)。q(t,x)為正表示從水到環境熱傳遞。

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

(8)

或:

(9)

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

θamb——環境溫度,可以是環境空氣溫度,也可以是地面溫度,℃

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

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

(10)

式中θout——控制體粒子在行程時間后到達出口時的當前溫度,℃

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

因此,管道中的控制體粒子溫度是漸近接近環境溫度的時間的指數函數,單個控制體粒子通過管道時的溫度分布見圖2。為了便于比較,采用線性插值法計算兩個相鄰斷點之間的溫度。當指數函數接近線性時,即當溫降較小時,線性近似是合理的。

圖2 單個控制體粒子通過管道時的溫度分布

采用漸變模式來表示時間區間內的入口溫度分布。即入口溫度在每個時間區間內逐漸變化。采用指數插值法時,入口溫度分布為分段指數函數,而采用線性插值法時,入口溫度分布為分段線性函數。

2.4 熱損失計算

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

在開始計算時間t0和結束計算時間tend之間,由環境造成的水的熱損失見式(11)～(14):

(11)

(12)

(13)

(14)

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

t0——開始計算時間,s

tend——結束計算時間,s

L——管道長度,m

Qin——入口能量,J

Qout——出口能量,J

U(tend)——tend時刻熱力學能,J

U(t0)——t0時刻熱力學能,J

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

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

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

U(t)——t時刻熱力學能,J

式(11)的右側項表示入口能量、出口能量和內部能量的變化。入口能量Qin和出口能量Qout是累積變量,在模擬過程中會進行更新。Qin和Qout可以分別更新。在產生新的水邊界時更新Qin,而在出口溫度函數的斷點處更新Qout。t時刻熱力學能U(t)是一個狀態變量,它只在t0和tend處計算。

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

在線性插值法中,假設入口溫度是t的分段線性函數,出口溫度分布用分段線性函數逼近。因為兩個斷點之間的流量恒定,所以積分在式(12)、(13)中由于線性被積函數的緣故可以被簡化。為了計算熱力學能,需要從入口到出口的水溫分布。作為位置x的函數的水溫分布也由分段線性函數來近似。

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

3 模型驗證

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

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

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

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

表1 案例供熱管道其他參數

① 出口溫度誤差

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

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

由圖4可以看出,除初始化階段外,出口溫度模擬值和測量值有較好的一致性。GL模型的最大誤差為±0.6 ℃。總體上看,出口溫度模擬值略低于測量值,這可能是對熱損失的輕微高估造成的,后續將進一步完善熱損失模型。這與特征線法和隱式迎風法[15]相當,優于一階隱式迎風法[19],但計算速度更快。

② 管道水溫分布

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

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

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

當水邊界對管道溫度分布影響很小時刪除水邊界是可能的,這涉及如何確定一個閾值,從而確定是否產生了新的事件。在大型模型中,刪除這種多余的影響不大的水邊界可能會顯著加快計算速度,刪除多余水邊界的標準是未來需要進一步研究的問題,這涉及計算量和精度的平衡。

4 結論及展望

① 提出了一種供熱管道動態熱力仿真方法。該方法基于拉格朗日方法,采用離散事件模擬(Discrete Event Simulation,DES)跟蹤沿管道傳輸的水力前沿,是首次將DES應用于供熱管道動態仿真的研究。該方法將水在管內的1個流動狀態變化視為1個事件,將整個供熱管網模擬分為多個連續事件的組合。事件可以被創建、修改、刪除。仿真的時間步長和空間步長可變,在達到相同計算精度的前提下,計算速度可大幅提升。

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

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

④ 通過研究同時發現,當水邊界對管道溫度分布影響很小時刪除水邊界是可能的,這涉及如何確定一個閾值,從而確定是否產生了新的事件。在大型模型中,刪除這種多余的影響不大的水邊界可能會顯著加快計算速度,刪除多余水邊界的標準是未來需要進一步研究的問題,這涉及計算量和精度的平衡。因此后續將繼續開發基于該核心思想的平臺產品,并研究該方法在大規模熱網和環狀熱網中的適用性和仿真速度。