基于Modelica与Julia的气—热—电综合能源系统动态特性分析

邱 勇，王 瑾，张树斌，何 院，张海明，丁 吉

（苏州同元软控信息技术有限公司，江苏 苏州 215000）

0 引言

综合能源系统在生产端整合了绿色清洁能源，通过输能网络与能源转换设备耦合，为消费端提供电、气、热/冷等形式的能源，可有效减少能源浪费，提升能源利用效率［1-2］。由于综合能源系统不同能流表现的流动特性不同，因此对其建模也与传统能源系统有一定差异。在综合能源系统模型中，既应包含同种能源传输网络的建模，也应包含异质能源之间相互转化、耦合的建模，增加了综合能源系统建模的困难程度［3-5］。

国内外学者对综合能源系统建模方法进行了大量研究。Geidl等［6］根据能源耦合元件的特性，明确提出了一种基于“ENERGY HUB”理论的综合能源系统集成建模方法，将系统内不同能量流耦合在一起，从而对综合能源系统运行中的耦合特性进行了描述。在此基础上，Moeini-Aghtaie 等［7］基于“ENERGY HUB”理论提出了综合能源系统的解耦方法，通过分析相关系数，从而详细介绍了综合能源系统的协同作用。Zhai 等［8］基于适用于双曲偏微分方程分析的特征方法，分析综合能源系统中天然气系统的动态变化，着重考虑了水锤压力引起的天然气回流问题，实现了其运行状态的真实模拟。Quelhas 等［9-10］提出了一种包括物理、经济、环境等方面的多周期广义网络流模型，以评估全系统能源流的经济效益。陈曦等［11-12］建立了综合考虑热量和天然气传输的动态特征的综合能源系统模型，该模型可用于综合能源系统仿真分析及不同能流子系统间的协同优化，为提升综合能源服务的经济性和安全性提供依据。Liu 等［13］对不同能源网络中多能互补的含义进行了定义，对包含电、气、热等3 种能源网络的综合能源系统进行综合分析，进而实现了对3 种能源系统的综合建模仿真。Wang等［14］在个体建模的细菌生态学建模方法中得到启发，认为在综合能源系统建模中应将大系统解耦为小个体。因此，他基于“个体为本”理论对综合能源系统进行建模，该方法将整个系统解耦为多个独立的个体，通过输入和输出集将每个个体统一交互。

综上，综合能源系统建模研究虽然已取得诸多进展，但对能流在网络系统中动态传输过程的模拟还相对较少，尤其是特殊流动现象还无法得到很好的模拟，而这在综合能源系统关键运行特性研究以及涉及到故障工况的建模研究中尤为重要。因此，为了研究综合能源系统运行时的动态特性，提出了基于Modelica 和Julia 语言的综合能源系统仿真分析方法。所提出的方法通过MWORKS 平台实现Modelica 语言与Julia 语言结合，将综合能源系统解耦为气、热、电3 个子系统，其中天然气、供热子系统通过Modelica 语言实现仿真，电力子系统通过Julia 语言实现仿真，有效避免了气、热、电等能流在网络中传输的时间尺度差异问题。最终将该方法应用于综合能源系统多种异质能流间的动态响应分析中，为综合能源系统运行优化和多主体在线能流分析提供理论支撑。

1 MWORKS平台与Modelica、Julia语言

1.1 MWORKS平台

如图1 所示，MWORKS 是基于国际知识统一表达与互联标准打造的系统智能设计与验证平台，包括科学计算环境-MWORKS.Syslab、系统建模仿真环境-MWORKS.Sysplorer、系统协同建模与模型数据管理平台-MWORKS.Syslink、工具箱-MWORKS.Toolbox 以及多领域工业模型库-MWORKS.Library 等。MWORKS 可以广泛应用于液压、传动、电气、热流、控制、动力学等多个专业领域，以及航天、航空、能源、车辆、工程机械等行业，支持系统/子系统/单机的知识积累、建模仿真、设计优化与运行维护需求［15-16］。

1.2 Modelica语言规范

Modelica 是一门建立在非因果建模思想上的、面向对象的、以数学方程为基础的物理系统建模语言，旨在支持连续和离散系统建模仿真计算。它采用基于广义基尔霍夫原理的连接机制进行统一建模，可以满足多领域需求。相较于其他计算机语言，Modelica 语言方便建模知识重用，为工程师提供便利。图2 为基于Modelica 的建模过程示意图，其主要有两种建模方式：①采用图形化工具建模，通过将Modelica 模型库的组件模型拖拽到建模视图区，根据实际结构关系绘制组件模型间的连接，然后给出各组件的参数，以建立整个系统的模型；②通过Modelica 语言编辑文本框的形式建模，在文本框内编辑各组件模型的代码以进行建模［17-19］。

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

1.3 Julia语言规范

Julia 是一个面向科学计算的高性能动态高级程序设计语言，其语法与其他科学计算语言相似。在许多情况下拥有能与编译型语言相媲美的性能。Julia 是个灵活的动态语言，适合科学和数值计算，性能可与传统静态类型语言媲美［20-21］。

Julia 语言便于与其他编程语言连接。由于Julia 还是一种相对较新的编程语言，这导致在特定领域缺少可用性的软件包。因此，Julia 语言具有一个与其他语言代码接口连接的能力，特别是使用其他编程语言提供的大量第三方库。同时，由其他编程语言编写的代码可以方便地转换成Julia 的代码，从而获得更具可读性和高效性的计算程序。如图3 所示，在MWORKS 平台，Julia 语言也可以与Modelica语言进行联合仿真，以达到协同建模的目的。

Fig.3 Joint simulation of Julia language and modelica language based on MWORKS platform图3 基于MWORKS平台的Julia语言与Modelica语言联合仿真示意图

2 综合能源系统建模方法

综合能源系统将气、热、电等多种能源系统耦合在一起协调运行，不同子系统间存在显著区别，因此在分析综合能源系统动态过程中应先建立各子系统的数学模型，进一步结合耦合模型，将各类设备模型和网络模型耦合在一起形成多能耦合的综合能源系统模型。

2.1 天然气系统模型

天然气在管道内的流动始终遵循质量守恒、动量守恒和能量守恒三大定律。考虑天然气管网自身特点，本文进行如下简化：①不考虑燃气在管道内流动过程与外界的换热；②管壁的膨胀可以忽略；③应用一元流动关系式；④管道截面温度均匀分布。

综上，气体流动基本方程可以写为：

连续性方程：

其中，A为管道截面积；ρ为燃气密度；w为燃气流速。

运动方程：

其中，λ为燃气管道的摩擦阻力系数；D为管道内径；α为管道倾斜角度；g为重力加速度。

能量方程：

其中，Q为气体的换热量；u为热力学能；H为焓。

为了对气体流动进行求解，采用交错网格形式的有限体积法［22］进行离散的数值方案。有限体积法的基本思路是：将计算区域划分为网格，并使每一个网格点周围有一个互不重复的控制体积；将待解微分方程（控制方程）对每一个控制体积积分，从而得出一组离散方程。交错网格形式是指将速度场和压力场的离散网格错开，不再重合在一起，如图4所示。

Fig.4 Schematic diagram of finite volume method using staggered grids图4 有限体积法—显示交错网格方法示意图

综上，基于交错网格的有限体积法管道流动模型可以写为：

质量守恒方程：

其中，m为单元格热媒的质量为单元格热媒的质量流量。

动量守恒方程：

其中，I为管道的冲量为管道的动量；Fg为重力的冲量；Ff为摩擦力的冲量；Fp为压差的冲量。

能量守恒方程：

式中，u为热水的比内能；h为比焓。

2.2 供热系统模型

为了研究供热系统管道内热量传递过程的动态变化，在供热管道水力计算模型基础上，耦合热力计算模型，形成供热管道水力—热力耦合的动态模型。供热管网热量传输变化可以表示为系统中热水的温度和质量流量的时空分布［23］，在不影响计算结果的情况下提出如下假设：①管道各部分均采用均质材料制成。热量仅沿管径传递，认为保温层和土层沿轴向无热传导；②管道与保温层、保温层与土壤视为无空隙接触，因而忽略了接触热阻的影响；③管道截面温度均匀分布；④忽略由水泵机械能耗散引起的热水温度上升。

图5 表示热水中的热量沿轴向传输和沿径向散失的过程。将管道离散为若干个单元，每个单元都有3 个关联的热流，这些热流决定了每个单元储存能量的变化。每个单元的径向传热都是独立的，通过管壁、保温层将热量散失到外界环境。轴向流动是整个管段径向温度场的连接条件，即每个单元的输入是上一个单元的输出，输出是下一个单元的输入。为了描述同时考虑轴向流动和径向散热的情况，建立了供热管道水力—热力耦合的动态模型。其中管道流动的连续性方程和运动方程与燃气管道相同，如式（4）和式（5）所示，能量方程在原来的基础上考虑管道的径向散热，如式（7）所示。

Fig.5 Physical model of pipe heat transfer图5 管道传热物理模型

其中，Ti，τ为热水温度；T0是环境温度；Rall表示热水和外部环境之间的总热阻。

2.3 电力系统潮流仿真模型

电力系统的基础组成部分主要包括：发电机、变压器、输电线路及负荷等。相较于燃气系统、供热系统，电力系统对于扰动的响应快速，一般在毫秒级，因此本文忽略了电力系统潮流的动态变化过程。在电力系统潮流仿真中，可根据节点电压方程建立电力系统潮流仿真模型，进而实现对电力系统各状态参数的求解［24］。

本文采用经典的电网潮流计算模型，其中潮流计算的基本方程为：

其中，Gij为节点i与节点j之间的电导；Bij为节点i与节点j之间的电纳；θij为节点i与节点j之间的相角。

2.4 耦合设备模型

热电联产（CHP）机组由燃气机组、余热锅炉、内燃机和外燃机等组成，是一种将供热和发电联合在一起的既产电又产热的生产方式［25］，相较于热电分产具有诸多优势，如降低能耗、提高空气质量、降低CO2排放、便于综合利用等。

CHP 机组其热功率与电功率之间的关系可表示为：

CHP 机组的耗气量为：

CHP 机组的向供热系统输出热量，其热出力与水温满足：

3 模型验证及模拟结果分析

3.1 模型验证

本文以气、热、电耦合的综合能源系统为研究对象，以气、热等能源传输中的动态特性为侧重点展开，因此在模型验证中主要验证燃气传输时的快瞬变流、慢瞬变流过程，水力传输过程的“水锤”现象，热力传输的“热惯性”等。

3.1.1 输气管道动态结果验证

对输气管道慢瞬变流工况进行测试［21］，给定出口流量，入口处压力随时间正弦变化，管段末端压力和管段输入流量变化如图6 所示，可以看到，模拟结果与实验结果在趋势上基本一致，在数值上误差较小，最大相对误差约为2.34%。

Fig.6 Verification of slow transient flow in gas pipelines图6 输气管道慢瞬变流验证

快瞬变流是指管道的突发状况，导致流动状态快速变化的情况。假设出口流量从1 kg/s 快速上升至788 kg/s，在30 min 时又快速下降至78.8 kg/s，仿真结果如图7 所示，模拟结果与实验结果在趋势上基本一致。由于实验值是在理想状态下得到的，模型中考虑了介质的物性，导致质量流率的结果误差较大。

Fig.7 Verification of fast transient flow in gas pipelines图7 输气管道快瞬变流验证

3.1.2 供热管道动态结果验证

图8 为供热管道质量流量突变时的动态水力特性变化曲线，以Modelica 标准库管道为对比对象，对供热管道进行测试。可以看到，模拟结果与对比管道结果基本一致，说明本文提出方法可以准确地模拟水力工况中的“水锤”现象。

Fig.8 Verification of dynamic hydraulic characteristics of heating pipelines图8 供热管道动态水力特性验证

供热管道热源温度突变时的动态热力特性变化曲线如图9 所示。可以看到，两者模拟结果吻合较好，最大相对误差为0.12%，产生误差的原因是模型中考虑了管道散热，导致温度略低。此外，由于热惯性的影响，各段的温度变化相对滞后。

Fig.9 Verification of dynamic thermodynamic characteristics of heating Pipelines图9 供热管道动态热力特性验证

3.2 综合能源系统动态特性分析

3.2.1 系统结构

本文构建了如图10 所示的综合能源系统进行案例研究，该综合能源系统主要包含电、气、热等能源子系统和CHP 等耦合设备。其中，电力系统为IEEE30 节点电力系统，通过CHP 机组与17 节点天然气系统和20 节点供热系统耦合在一起。

Fig.10 Topology of integrated energy system图10 综合能源系统拓扑结构

在上述综合能源系统模型和综合能源系统拓扑结构的基础上，基于Modelica 和Julia 语言完成对上述综合能源系统组件模型开发、系统模型搭建及调试，综合能源系统模型如图11 所示。其中，气网和热网的部分基于Modelica语言完成，电网部分基于Julia 语言完成，并在MWORKS 平台下完成联合仿真，实现电-气-热综合能源系统统一求解。

Fig.11 Integrated energy system model图11 综合能源系统模型

3.2.2 负荷变化工况仿真结果及分析

由于某些地区环境温度波动较大，需要调节供热量以保证用户室温处于适宜水平。为探究这种由于热负荷变化引起系统的状态变化情况，假设供热系统在运行2 min后供热量逐渐下降10%，其他条件不变，其中CHP 机组处于“以热定电”模式运行。

图12 为供热系统S1 管网各供水管出口温度随时间的变化曲线。可以看到，供热量逐渐下降10%后各供水管出口温度也相继下降，由于热惯性的影响，下降时间约为30～50 min，最终稳定温度约为80.9 ℃、80.8 ℃、80.7 ℃、80.6 ℃、80.4 ℃。同时，由于水力传输与热力传输的耦合，从第1 根管道响应到第5 根管道响应，两者之间的时间差约为30 min，这体现了热量在管道内传输的延时特性。在综合能源系统的优化调度中，对于热传输延时的描述至关重要，热传输延时不仅直接影响系统的热能传输效率和响应速度，而且会严重影响供热量与热负荷之间的供需平衡关系。

Fig.12 Temperature changes of the inlet and outlet of each water supply pipe in the heating system图12 供热系统各供水管出口温度变化

由于CHP 机组在“以热定电”模式运行，当供热量下降10%时，其热、电出力也相应降低，此时电力系统电负荷需求由平衡节点承担，平衡节点的电功率变化和耗气量变化如图13 所示。可以看到，为满足电负荷需求，平衡节点电压由0.52 pu 上升至0.71 pu，耗气量由2.55 kg/s 上升至3.51 kg/s。电力系统和供热系统之间的相互耦合，导致供热系统内部扰动作用于电力系统的源侧，引起了电能输出的重新分配。

Fig.13 Changes of electric power and gas consumption at balancing nodes in power systems图13 电力系统平衡节点电功率及其耗气量变化

电力系统和供热系统的状态变化也同时引起了天然气系统的响应，首先在燃气系统的荷侧产生影响，导致用气负荷快速变化，进而引起各支管的压力和传输流量的重新分配。燃气系统各管段出口的压力和流量变化如图14和图15 所示，可以看到，各管段末端压力和出口流量均会受影响而产生波动，并在一段时间（大约10～20 min）后再次达到稳定状态。由于天然气系统连接的各终端用户不同，不同区域管段末端压力和出口流量的变化也不同，其中，管段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 燃气系统各管段末端压力变化

Fig.15 Flow changes of the inlet and outlet of each water supply pipe in the heating system图15 燃气系统各管段出口流量变化

上述研究表明，通过基于Modelica 与Julia 语言的综合能源系统仿真分析方法，既能够模拟单一能源的网络动态特性，也能够模拟多种异质能流间的动态响应，且模拟效果良好，该方法可为综合能源系统运行优化和多主体在线能流分析等研究提供理论支持。

4 研究结论

本文提出了考虑网络动态特性的综合能源系统数学模型，并以此为基础，通过对综合能源系统变负荷工况的仿真，深入探讨了综合能源系统运行的动态特性，主要结论如下：

（1）基于综合能源系统模型对气、热等能源传输中的动态特性进行模拟验证，结果表明，本文所提出方法对模拟输气管道的快瞬变流、慢瞬变流特性以及传输过程的“水锤”“热惯性”等现象具有良好效果。其中，对输气管道慢瞬变流特性的仿真误差小于2.34%，对“水锤”“热惯性”等现象的仿真误差小于0.12%，验证了方法的可靠性。

（2）由于综合能源系统的多能耦合特性，供热系统供热量下降10%后，首先会引起供热网络的状态变化，并在30～50 min 后再次达到稳定状态。同时，供热量的变化会通过耦合单元引起耦合子系统的响应，进而引起整个综合能源系统潮流重新分布。相较于供热系统，电力系统和天然气系统响应较快，其中天然气系统的响应时间约为10～20 min。

（3）本研究可应用于综合能源系统运行优化与多主体在线能流分析，为提高综合能源系统的科学管理能力提供理论依据。