邱 勇，王 瑾，張樹斌，何 院，張海明，丁 吉

（蘇州同元軟控信息技術有限公司，江蘇 蘇州 215000）

0 引言

綜合能源系統(tǒng)在生產端整合了綠色清潔能源，通過輸能網絡與能源轉換設備耦合，為消費端提供電、氣、熱/冷等形式的能源，可有效減少能源浪費，提升能源利用效率［1-2］。由于綜合能源系統(tǒng)不同能流表現(xiàn)的流動特性不同，因此對其建模也與傳統(tǒng)能源系統(tǒng)有一定差異。在綜合能源系統(tǒng)模型中，既應包含同種能源傳輸網絡的建模，也應包含異質能源之間相互轉化、耦合的建模，增加了綜合能源系統(tǒng)建模的困難程度［3-5］。

國內外學者對綜合能源系統(tǒng)建模方法進行了大量研究。Geidl等［6］根據能源耦合元件的特性，明確提出了一種基于“ENERGY HUB”理論的綜合能源系統(tǒng)集成建模方法，將系統(tǒng)內不同能量流耦合在一起，從而對綜合能源系統(tǒng)運行中的耦合特性進行了描述。在此基礎上，Moeini-Aghtaie 等［7］基于“ENERGY HUB”理論提出了綜合能源系統(tǒng)的解耦方法，通過分析相關系數(shù)，從而詳細介紹了綜合能源系統(tǒng)的協(xié)同作用。Zhai 等［8］基于適用于雙曲偏微分方程分析的特征方法，分析綜合能源系統(tǒng)中天然氣系統(tǒng)的動態(tài)變化，著重考慮了水錘壓力引起的天然氣回流問題，實現(xiàn)了其運行狀態(tài)的真實模擬。Quelhas 等［9-10］提出了一種包括物理、經濟、環(huán)境等方面的多周期廣義網絡流模型，以評估全系統(tǒng)能源流的經濟效益。陳曦等［11-12］建立了綜合考慮熱量和天然氣傳輸?shù)膭討B(tài)特征的綜合能源系統(tǒng)模型，該模型可用于綜合能源系統(tǒng)仿真分析及不同能流子系統(tǒng)間的協(xié)同優(yōu)化，為提升綜合能源服務的經濟性和安全性提供依據。Liu 等［13］對不同能源網絡中多能互補的含義進行了定義，對包含電、氣、熱等3 種能源網絡的綜合能源系統(tǒng)進行綜合分析，進而實現(xiàn)了對3 種能源系統(tǒng)的綜合建模仿真。Wang等［14］在個體建模的細菌生態(tài)學建模方法中得到啟發(fā)，認為在綜合能源系統(tǒng)建模中應將大系統(tǒng)解耦為小個體。因此，他基于“個體為本”理論對綜合能源系統(tǒng)進行建模，該方法將整個系統(tǒng)解耦為多個獨立的個體，通過輸入和輸出集將每個個體統(tǒng)一交互。

綜上，綜合能源系統(tǒng)建模研究雖然已取得諸多進展，但對能流在網絡系統(tǒng)中動態(tài)傳輸過程的模擬還相對較少，尤其是特殊流動現(xiàn)象還無法得到很好的模擬，而這在綜合能源系統(tǒng)關鍵運行特性研究以及涉及到故障工況的建模研究中尤為重要。因此，為了研究綜合能源系統(tǒng)運行時的動態(tài)特性，提出了基于Modelica 和Julia 語言的綜合能源系統(tǒng)仿真分析方法。所提出的方法通過MWORKS 平臺實現(xiàn)Modelica 語言與Julia 語言結合，將綜合能源系統(tǒng)解耦為氣、熱、電3 個子系統(tǒng)，其中天然氣、供熱子系統(tǒng)通過Modelica 語言實現(xiàn)仿真，電力子系統(tǒng)通過Julia 語言實現(xiàn)仿真，有效避免了氣、熱、電等能流在網絡中傳輸?shù)臅r間尺度差異問題。最終將該方法應用于綜合能源系統(tǒng)多種異質能流間的動態(tài)響應分析中，為綜合能源系統(tǒng)運行優(yōu)化和多主體在線能流分析提供理論支撐。

1 MWORKS平臺與Modelica、Julia語言

1.1 MWORKS平臺

如圖1 所示，MWORKS 是基于國際知識統(tǒng)一表達與互聯(lián)標準打造的系統(tǒng)智能設計與驗證平臺，包括科學計算環(huán)境-MWORKS.Syslab、系統(tǒng)建模仿真環(huán)境-MWORKS.Sysplorer、系統(tǒng)協(xié)同建模與模型數(shù)據管理平臺-MWORKS.Syslink、工具箱-MWORKS.Toolbox 以及多領域工業(yè)模型庫-MWORKS.Library 等。MWORKS 可以廣泛應用于液壓、傳動、電氣、熱流、控制、動力學等多個專業(yè)領域，以及航天、航空、能源、車輛、工程機械等行業(yè)，支持系統(tǒng)/子系統(tǒng)/單機的知識積累、建模仿真、設計優(yōu)化與運行維護需求［15-16］。

1.2 Modelica語言規(guī)范

Modelica 是一門建立在非因果建模思想上的、面向對象的、以數(shù)學方程為基礎的物理系統(tǒng)建模語言，旨在支持連續(xù)和離散系統(tǒng)建模仿真計算。它采用基于廣義基爾霍夫原理的連接機制進行統(tǒng)一建模，可以滿足多領域需求。相較于其他計算機語言，Modelica 語言方便建模知識重用，為工程師提供便利。圖2 為基于Modelica 的建模過程示意圖，其主要有兩種建模方式：①采用圖形化工具建模，通過將Modelica 模型庫的組件模型拖拽到建模視圖區(qū)，根據實際結構關系繪制組件模型間的連接，然后給出各組件的參數(shù)，以建立整個系統(tǒng)的模型；②通過Modelica 語言編輯文本框的形式建模，在文本框內編輯各組件模型的代碼以進行建模［17-19］。

Fig.2 Schematic diagram of modeling process based on Modelica圖2 基于Modelica的建模過程示意圖

1.3 Julia語言規(guī)范

Julia 是一個面向科學計算的高性能動態(tài)高級程序設計語言，其語法與其他科學計算語言相似。在許多情況下?lián)碛心芘c編譯型語言相媲美的性能。Julia 是個靈活的動態(tài)語言，適合科學和數(shù)值計算，性能可與傳統(tǒng)靜態(tài)類型語言媲美［20-21］。

Julia 語言便于與其他編程語言連接。由于Julia 還是一種相對較新的編程語言，這導致在特定領域缺少可用性的軟件包。因此，Julia 語言具有一個與其他語言代碼接口連接的能力，特別是使用其他編程語言提供的大量第三方庫。同時，由其他編程語言編寫的代碼可以方便地轉換成Julia 的代碼，從而獲得更具可讀性和高效性的計算程序。如圖3 所示，在MWORKS 平臺，Julia 語言也可以與Modelica語言進行聯(lián)合仿真，以達到協(xié)同建模的目的。

Fig.3 Joint simulation of Julia language and modelica language based on MWORKS platform圖3 基于MWORKS平臺的Julia語言與Modelica語言聯(lián)合仿真示意圖

2 綜合能源系統(tǒng)建模方法

綜合能源系統(tǒng)將氣、熱、電等多種能源系統(tǒng)耦合在一起協(xié)調運行，不同子系統(tǒng)間存在顯著區(qū)別，因此在分析綜合能源系統(tǒng)動態(tài)過程中應先建立各子系統(tǒng)的數(shù)學模型，進一步結合耦合模型，將各類設備模型和網絡模型耦合在一起形成多能耦合的綜合能源系統(tǒng)模型。

2.1 天然氣系統(tǒng)模型

天然氣在管道內的流動始終遵循質量守恒、動量守恒和能量守恒三大定律?？紤]天然氣管網自身特點，本文進行如下簡化：①不考慮燃氣在管道內流動過程與外界的換熱；②管壁的膨脹可以忽略；③應用一元流動關系式；④管道截面溫度均勻分布。

綜上，氣體流動基本方程可以寫為：

連續(xù)性方程：

其中，A為管道截面積；ρ為燃氣密度；w為燃氣流速。

運動方程：

其中，λ為燃氣管道的摩擦阻力系數(shù)；D為管道內徑；α為管道傾斜角度；g為重力加速度。

能量方程：

其中，Q為氣體的換熱量；u為熱力學能；H為焓。

為了對氣體流動進行求解，采用交錯網格形式的有限體積法［22］進行離散的數(shù)值方案。有限體積法的基本思路是：將計算區(qū)域劃分為網格，并使每一個網格點周圍有一個互不重復的控制體積；將待解微分方程（控制方程）對每一個控制體積積分，從而得出一組離散方程。交錯網格形式是指將速度場和壓力場的離散網格錯開，不再重合在一起，如圖4所示。

Fig.4 Schematic diagram of finite volume method using staggered grids圖4 有限體積法—顯示交錯網格方法示意圖

綜上，基于交錯網格的有限體積法管道流動模型可以寫為：

質量守恒方程：

其中，m為單元格熱媒的質量為單元格熱媒的質量流量。

動量守恒方程：

其中，I為管道的沖量為管道的動量；Fg為重力的沖量；Ff為摩擦力的沖量；Fp為壓差的沖量。

能量守恒方程：

式中，u為熱水的比內能；h為比焓。

2.2 供熱系統(tǒng)模型

為了研究供熱系統(tǒng)管道內熱量傳遞過程的動態(tài)變化，在供熱管道水力計算模型基礎上，耦合熱力計算模型，形成供熱管道水力—熱力耦合的動態(tài)模型。供熱管網熱量傳輸變化可以表示為系統(tǒng)中熱水的溫度和質量流量的時空分布［23］，在不影響計算結果的情況下提出如下假設：①管道各部分均采用均質材料制成。熱量僅沿管徑傳遞，認為保溫層和土層沿軸向無熱傳導；②管道與保溫層、保溫層與土壤視為無空隙接觸，因而忽略了接觸熱阻的影響；③管道截面溫度均勻分布；④忽略由水泵機械能耗散引起的熱水溫度上升。

圖5 表示熱水中的熱量沿軸向傳輸和沿徑向散失的過程。將管道離散為若干個單元，每個單元都有3 個關聯(lián)的熱流，這些熱流決定了每個單元儲存能量的變化。每個單元的徑向傳熱都是獨立的，通過管壁、保溫層將熱量散失到外界環(huán)境。軸向流動是整個管段徑向溫度場的連接條件，即每個單元的輸入是上一個單元的輸出，輸出是下一個單元的輸入。為了描述同時考慮軸向流動和徑向散熱的情況，建立了供熱管道水力—熱力耦合的動態(tài)模型。其中管道流動的連續(xù)性方程和運動方程與燃氣管道相同，如式（4）和式（5）所示，能量方程在原來的基礎上考慮管道的徑向散熱，如式（7）所示。

Fig.5 Physical model of pipe heat transfer圖5 管道傳熱物理模型

其中，Ti，τ為熱水溫度；T0是環(huán)境溫度；Rall表示熱水和外部環(huán)境之間的總熱阻。

2.3 電力系統(tǒng)潮流仿真模型

電力系統(tǒng)的基礎組成部分主要包括：發(fā)電機、變壓器、輸電線路及負荷等。相較于燃氣系統(tǒng)、供熱系統(tǒng)，電力系統(tǒng)對于擾動的響應快速，一般在毫秒級，因此本文忽略了電力系統(tǒng)潮流的動態(tài)變化過程。在電力系統(tǒng)潮流仿真中，可根據節(jié)點電壓方程建立電力系統(tǒng)潮流仿真模型，進而實現(xiàn)對電力系統(tǒng)各狀態(tài)參數(shù)的求解［24］。

本文采用經典的電網潮流計算模型，其中潮流計算的基本方程為：

其中，Gij為節(jié)點i與節(jié)點j之間的電導；Bij為節(jié)點i與節(jié)點j之間的電納；θij為節(jié)點i與節(jié)點j之間的相角。

2.4 耦合設備模型

熱電聯(lián)產（CHP）機組由燃氣機組、余熱鍋爐、內燃機和外燃機等組成，是一種將供熱和發(fā)電聯(lián)合在一起的既產電又產熱的生產方式［25］，相較于熱電分產具有諸多優(yōu)勢，如降低能耗、提高空氣質量、降低CO2排放、便于綜合利用等。

CHP 機組其熱功率與電功率之間的關系可表示為：

CHP 機組的耗氣量為：

CHP 機組的向供熱系統(tǒng)輸出熱量，其熱出力與水溫滿足：

3 模型驗證及模擬結果分析

3.1 模型驗證

本文以氣、熱、電耦合的綜合能源系統(tǒng)為研究對象，以氣、熱等能源傳輸中的動態(tài)特性為側重點展開，因此在模型驗證中主要驗證燃氣傳輸時的快瞬變流、慢瞬變流過程，水力傳輸過程的“水錘”現(xiàn)象，熱力傳輸?shù)摹盁釕T性”等。

3.1.1 輸氣管道動態(tài)結果驗證

對輸氣管道慢瞬變流工況進行測試［21］，給定出口流量，入口處壓力隨時間正弦變化，管段末端壓力和管段輸入流量變化如圖6 所示，可以看到，模擬結果與實驗結果在趨勢上基本一致，在數(shù)值上誤差較小，最大相對誤差約為2.34%。

Fig.6 Verification of slow transient flow in gas pipelines圖6 輸氣管道慢瞬變流驗證

快瞬變流是指管道的突發(fā)狀況，導致流動狀態(tài)快速變化的情況。假設出口流量從1 kg/s 快速上升至788 kg/s，在30 min 時又快速下降至78.8 kg/s，仿真結果如圖7 所示，模擬結果與實驗結果在趨勢上基本一致。由于實驗值是在理想狀態(tài)下得到的，模型中考慮了介質的物性，導致質量流率的結果誤差較大。

Fig.7 Verification of fast transient flow in gas pipelines圖7 輸氣管道快瞬變流驗證

3.1.2 供熱管道動態(tài)結果驗證

圖8 為供熱管道質量流量突變時的動態(tài)水力特性變化曲線，以Modelica 標準庫管道為對比對象，對供熱管道進行測試?？梢钥吹剑M結果與對比管道結果基本一致，說明本文提出方法可以準確地模擬水力工況中的“水錘”現(xiàn)象。

Fig.8 Verification of dynamic hydraulic characteristics of heating pipelines圖8 供熱管道動態(tài)水力特性驗證

供熱管道熱源溫度突變時的動態(tài)熱力特性變化曲線如圖9 所示?？梢钥吹?，兩者模擬結果吻合較好，最大相對誤差為0.12%，產生誤差的原因是模型中考慮了管道散熱，導致溫度略低。此外，由于熱慣性的影響，各段的溫度變化相對滯后。

Fig.9 Verification of dynamic thermodynamic characteristics of heating Pipelines圖9 供熱管道動態(tài)熱力特性驗證

3.2 綜合能源系統(tǒng)動態(tài)特性分析

3.2.1 系統(tǒng)結構

本文構建了如圖10 所示的綜合能源系統(tǒng)進行案例研究，該綜合能源系統(tǒng)主要包含電、氣、熱等能源子系統(tǒng)和CHP 等耦合設備。其中，電力系統(tǒng)為IEEE30 節(jié)點電力系統(tǒng)，通過CHP 機組與17 節(jié)點天然氣系統(tǒng)和20 節(jié)點供熱系統(tǒng)耦合在一起。

Fig.10 Topology of integrated energy system圖10 綜合能源系統(tǒng)拓撲結構

在上述綜合能源系統(tǒng)模型和綜合能源系統(tǒng)拓撲結構的基礎上，基于Modelica 和Julia 語言完成對上述綜合能源系統(tǒng)組件模型開發(fā)、系統(tǒng)模型搭建及調試，綜合能源系統(tǒng)模型如圖11 所示。其中，氣網和熱網的部分基于Modelica語言完成，電網部分基于Julia 語言完成，并在MWORKS 平臺下完成聯(lián)合仿真，實現(xiàn)電-氣-熱綜合能源系統(tǒng)統(tǒng)一求解。

Fig.11 Integrated energy system model圖11 綜合能源系統(tǒng)模型

3.2.2 負荷變化工況仿真結果及分析

由于某些地區(qū)環(huán)境溫度波動較大，需要調節(jié)供熱量以保證用戶室溫處于適宜水平。為探究這種由于熱負荷變化引起系統(tǒng)的狀態(tài)變化情況，假設供熱系統(tǒng)在運行2 min后供熱量逐漸下降10%，其他條件不變，其中CHP 機組處于“以熱定電”模式運行。

圖12 為供熱系統(tǒng)S1 管網各供水管出口溫度隨時間的變化曲線?？梢钥吹剑崃恐饾u下降10%后各供水管出口溫度也相繼下降，由于熱慣性的影響，下降時間約為30～50 min，最終穩(wěn)定溫度約為80.9 ℃、80.8 ℃、80.7 ℃、80.6 ℃、80.4 ℃。同時，由于水力傳輸與熱力傳輸?shù)鸟詈希瑥牡? 根管道響應到第5 根管道響應，兩者之間的時間差約為30 min，這體現(xiàn)了熱量在管道內傳輸?shù)难訒r特性。在綜合能源系統(tǒng)的優(yōu)化調度中，對于熱傳輸延時的描述至關重要，熱傳輸延時不僅直接影響系統(tǒng)的熱能傳輸效率和響應速度，而且會嚴重影響供熱量與熱負荷之間的供需平衡關系。

Fig.12 Temperature changes of the inlet and outlet of each water supply pipe in the heating system圖12 供熱系統(tǒng)各供水管出口溫度變化

由于CHP 機組在“以熱定電”模式運行，當供熱量下降10%時，其熱、電出力也相應降低，此時電力系統(tǒng)電負荷需求由平衡節(jié)點承擔，平衡節(jié)點的電功率變化和耗氣量變化如圖13 所示。可以看到，為滿足電負荷需求，平衡節(jié)點電壓由0.52 pu 上升至0.71 pu，耗氣量由2.55 kg/s 上升至3.51 kg/s。電力系統(tǒng)和供熱系統(tǒng)之間的相互耦合，導致供熱系統(tǒng)內部擾動作用于電力系統(tǒng)的源側，引起了電能輸出的重新分配。

Fig.13 Changes of electric power and gas consumption at balancing nodes in power systems圖13 電力系統(tǒng)平衡節(jié)點電功率及其耗氣量變化

電力系統(tǒng)和供熱系統(tǒng)的狀態(tài)變化也同時引起了天然氣系統(tǒng)的響應，首先在燃氣系統(tǒng)的荷側產生影響，導致用氣負荷快速變化，進而引起各支管的壓力和傳輸流量的重新分配。燃氣系統(tǒng)各管段出口的壓力和流量變化如圖14和圖15 所示，可以看到，各管段末端壓力和出口流量均會受影響而產生波動，并在一段時間（大約10～20 min）后再次達到穩(wěn)定狀態(tài)。由于天然氣系統(tǒng)連接的各終端用戶不同，不同區(qū)域管段末端壓力和出口流量的變化也不同，其中，管段2-5 的管段末端壓力上升，出口流量呈下降趨勢；管段6-8 的管段末端壓力下降，出口流量呈上升趨勢；管段9-16的管段末端壓力下降，出口流量短暫波動后恢復。

Fig.14 Pressure changes of the inlet and outlet of each water supply pipe in the heating system圖14 燃氣系統(tǒng)各管段末端壓力變化

Fig.15 Flow changes of the inlet and outlet of each water supply pipe in the heating system圖15 燃氣系統(tǒng)各管段出口流量變化

上述研究表明，通過基于Modelica 與Julia 語言的綜合能源系統(tǒng)仿真分析方法，既能夠模擬單一能源的網絡動態(tài)特性，也能夠模擬多種異質能流間的動態(tài)響應，且模擬效果良好，該方法可為綜合能源系統(tǒng)運行優(yōu)化和多主體在線能流分析等研究提供理論支持。

4 研究結論

本文提出了考慮網絡動態(tài)特性的綜合能源系統(tǒng)數(shù)學模型，并以此為基礎，通過對綜合能源系統(tǒng)變負荷工況的仿真，深入探討了綜合能源系統(tǒng)運行的動態(tài)特性，主要結論如下：

（1）基于綜合能源系統(tǒng)模型對氣、熱等能源傳輸中的動態(tài)特性進行模擬驗證，結果表明，本文所提出方法對模擬輸氣管道的快瞬變流、慢瞬變流特性以及傳輸過程的“水錘”“熱慣性”等現(xiàn)象具有良好效果。其中，對輸氣管道慢瞬變流特性的仿真誤差小于2.34%，對“水錘”“熱慣性”等現(xiàn)象的仿真誤差小于0.12%，驗證了方法的可靠性。

（2）由于綜合能源系統(tǒng)的多能耦合特性，供熱系統(tǒng)供熱量下降10%后，首先會引起供熱網絡的狀態(tài)變化，并在30～50 min 后再次達到穩(wěn)定狀態(tài)。同時，供熱量的變化會通過耦合單元引起耦合子系統(tǒng)的響應，進而引起整個綜合能源系統(tǒng)潮流重新分布。相較于供熱系統(tǒng)，電力系統(tǒng)和天然氣系統(tǒng)響應較快，其中天然氣系統(tǒng)的響應時間約為10～20 min。

（3）本研究可應用于綜合能源系統(tǒng)運行優(yōu)化與多主體在線能流分析，為提高綜合能源系統(tǒng)的科學管理能力提供理論依據。