Dynamic simulation of GEH-IES with distributed parameter characteristics for hydrogen-blending transportation

Abstract

For the purpose of environment protecting and energy saving, renewable energy has been distributed into the power grid in a considerable scale. However, the consuming capacity of the power grid for renewable energy is relatively limited. As an effective way to absorb the excessive renewable energy, the power to gas (P2G) technology is able to convert excessive renewable energy into hydrogen. Hydrogen-blending natural gas pipeline is an efficient approach for hydrogen transportation. However, hydrogen-blending natural gas complicates the whole integrated energy system (IES), making it more problematic to cope with the equipment failure, demand response and dynamic optimization. Nevertheless, dynamic simulation of distribution parameters of gas–electricity–hydrogen (GEH) energy system, especially for hydrogen concentration, still remains a challenge. The dynamics of hydrogen-blending IES is undiscovered. To tackle the issue, an iterative solving framework of the GEH-IES and a cell segment-based method for hydrogen mixing ratio distribution are proposed in this paper. Two typical numerical cases studying the conditions under which renewables fluctuate and generators fail are conducted on a real-word system. The results show that hydrogen blending timely and spatially influences the flow parameters, of which the hydrogen mixing ratio and gas pressure loss along the gas pipeline are negatively correlated and the response to hydrogen mixing ratio is time-delayed. Moreover, the hydrogen-blending amount and position also have a significant impact on the performance of the compressor.

Graphical abstract

Keywords

gas–electricity IES / dynamic simulation / hydrogen blending / power to gas (P2G) / renewable energy

1 Introduction

For the purpose of protecting environment, the proportion of primary energy all over the world is gradually changing with use of more renewable energy and less traditional fossil energy [1]. However, renewable energy, mostly excessive for the system, is discontinuous and instable, making the system difficult to match the user load timely. Integrated energy system (IES) with coordinated operation of various energy sources and distribution networks is effective to overcome such discontinuity and instability [2,3], and able to timely consume the excessive renewable energy and deliver the stored energy when users need [4]. Power to gas (P2G) [5] is a typical technology to realize it in IESs with a considerable storage density. Compared with traditional energy, hydrogen energy has clean, sustainable advantages, and can replace traditional fossil energy in many fields, i.e., aerospace flight engine fuel [6]. However, constructing the needed gas storage facility and transporting hydrogen are cost-expensive, safety-problematic, and technology-challenging [7]. Recently, as a new rising technology, hydrogen-blending natural gas pipeline transportation shows a promising potential to achieve the lower cost, lower risk, and less technical obstacles relatively with far less storage and transportation facilities [8,9].

The gas–electricity–hydrogen IES (GEH-IES) is an IES with a great potential that leverages the interconversion of multiple energy sources, including natural gas, electricity and hydrogen, to achieve more efficient, sustainable, and environmentally friendly energy production and utilization. To further guarantee the performance of operational economy and system safety, understanding of the dynamics of such a GEH-IES is of great significance [10]. Since the distribution of hydrogen concentration along the natural gas pipelines has a significant impact on the performance of equipment, solutions to distribution flow parameters and gas components are needed [11]. Along with the integration of the power grid, the natural gas network, and renewable energy sources, dynamic simulation of distribution parameters for the GEH energy system faces the challenge of high dimensionality, high instability, and complex coupling characteristics.

Actually, dynamic simulation and optimization of the gas–electricity IES (GE-IES) has been discussed for many years, and the study mainly focuses on multiple-time-scale issues for tackling the huge differences on the time scale between power grid and gas networks [12]. Shen et al. [13 ] proposed the method of singular perturbation to decompose GE-IES into fast and slow subsystems and obtain the solutions in three-time scales. From a perspective of decoupling GE-IES, Ma et al. [14] presented an optimization method to find an optimal scheduling plan. On the other hand, transient simulation of natural gas network with maximum inertia in GE-IES is studied with more attention. Numeric simulation methods are mainly focused on these issues. Since hydrogen blending brings new issues to the dynamic gas mixture process, along with the original thermodynamic, hydraulic and physical characteristics calculation, mathematical modeling of the natural gas pipeline system with hydrogen blending poses challenges. To have a deeper sight into the dynamics of the hydrogen-blending natural gas pipeline system, some modeling and simulation are also conducted. Tabkhi et al. [15] proposed a method for modeling and estimating the parameters of the gas pipeline system with hydrogen blending. Agaie et al. [16] proposed the reduced-order model for solving transient hydrogen–natural gas flow. Guandalini et al. [17,18] presented a transient model and simulated the quality of the hydrogen-blending natural gas along the pipeline. Considering the load of industrial and residential users, Cheli et al. [19] built a simulation model for a low-pressure natural gas distribution network with single local hydrogen injection point while Elaoud & Hadj-Taïeb [20] conducted research on transient simulation of high-pressure hydrogen–natural gas flow along rigid pipelines. Liu et al. [21] performed an economic analysis for large-scale renewable hydrogen application on the basis of modeling natural gas pipelines. There are also several projects focusing on the engineering implementation of blending hydrogen in natural gas networks [22–24].

Although improvement on dynamic simulation of hydrogen–blending natural gas networks has been made to some extent, study on the overall solution scheme to tackle the hydrogen issues on GEH-IES is still in its infancy. Some researches focusing on the topic of modeling, simulation, and optimization of IES with hydrogen storage and utilization have been conducted. Colbertaldo et al. [25] built models of the IES with the P2G technology in Italy and analyzed the influence on transporting economy. Fang [26] conducted the life cycle estimation of wind power–hydrogen IES, proving its huge cut on the energy cost. Likewise, a considerable number of researches on the optimization planning of the IES with hydrogen storage and utilization are conducted, such as the wind–photovoltaic system [27–29], combined cooling, heating and power system (CCHP) [30] and GE-IES [31,32].

However, the aforementioned researches focus on the modeling and optimization of GEH-IES of hydrogen-storage of systems, with no consideration on the transportation of hydrogen via gas pipeline. Studies on hydrogen blending transportation along the gas pipeline in GEH-IES remains rare. Therefore, the dynamic transporting characteristics of the GEH-IES and its dynamics are not clearly discovered, although some scholars are starting to study this issue. Zhou et al. [33] conducted modeling and simulation of the GEH-IES and analyzed system characteristics with different hydrogen-injecting modes in static scenarios. In addition, Liu et al. [34] disclosed the influence of power-to-hydrogen equipment on the environment and economy based on static GE system modeling.

Moreover, researches on the modeling and analysis of GEH-IES also remain problematic. The model for GEH-IES is simplified and no transient simulation is considered, making the numerical cases far different from those of the real world. Simulation of distribution parameters, especially for hydrogen density and gas property, is not conducted, which plays a key role in detecting hydrogen embrittlement in gas pipelines [35–38]. Due to time-spatial evolution characteristics of hydrogen concentration, the instability of dynamics of natural gas network increases [39], leaving it unclear when faced with the fluctuation of renewable energy, equipment failures and users’ demands. The key to solving the above problems is to propose a transient simulation method to fulfill the distributed parameter solving in GEH-IES, mainly the hydrogen mixing rate and flow parameters, with consideration of the instable renewable energy power. In other words, a fine-grained time-space gas flow and gas property (or gas component) solving scheme cooperating with the traditional power flow optimization framework is urgently needed. However, there is still no literature focusing this topic (transient simulation of GEH-IES), the overall dynamic characteristics of GEH-IES remain undiscovered. Moreover, due to the different time-space resolution of the hydrogen-blending natural gas network and power gird, parameters on coupled nodes shall be mainly focused on in such transient simulation methods.

To tackle the above issues of GEH-IES, a novel iterative-interactive simulation scheme based on the cell-segment flow modeling method is proposed, making the simulation possible on the dynamic operating conditions and space-distribution parameters including gas flow mass, gas pressure, and gas property. Cooperating with power flow optimization, the dynamic evolution of GEH-IES parameters can be observed. Additionally, two typical transient cases are conducted to understand the dynamic behaviors of the system.

This paper is contributive because a novel iterative-interactive calculation scheme based on cell segment solution method with fine-grained resolution on time and space is presented and validated, providing a general framework for dynamic simulation and analysis of GEH-IES. In addition, the impact of hydrogen blending operation on flow parameters, hydrogen mixing ratio, and compressor performance is disclosed and responding dynamic characteristics of the gas networks is discovered.

2 Background

2.1 Modeling of natural gas network

In actual operation, natural gas is transported by pipeline networks composed of multiple gas transmission pipelines [40], gas supply stations, and gas compressor stations including electricity-driven or gas-driven compressors, where gas supply stations are represented as nodes and the gas transmission pipelines are edges [41,42].

The gas network is dominated by flow continuity equation and pressure balance equation of the flow loop [43].

\mathit{C}\mathit{q}=\mathit{Q},

(1)\mathit{D}\times \mathrm{\Delta }\mathit{p}=\mathbf{0},

where \mathit{C} is the node-edge adjacent matrix of pipe networks, \mathit{q} is the flow vector of edges including the flow of inlet and outlet nodes, \mathit{Q} is the load vector of nodes, \mathit{D} is loop-edge adjacent matrix, and \mathrm{\Delta }\mathit{p} is the drop vector of the pressure of edges [43].

The elements of \mathrm{\Delta }\mathit{p} are determined by the Weymouth equation [44], for static flow

q^2=e\sqrt{\left|p_0^2-p_{\mathrm{L}}^2\right|},

(2)\mathrm{\Delta }p=p_{\mathrm{L}}^{}-p_0^{},

where q is the mass flow of the pipe, p_0^{} and p_{\mathrm{L}}^{} are the inlet and outlet pressure of the pipe, e is a constant related to the property of gas and pipeline. Notice the subscripts 0 and L mean the inlet and outlet of the gas pipeline, respectively.

2.1.1 Model of natural gas pipe

Assuming the temperature and altitude (where the angle between the pipeline and the horizontal plane is 0) along the pipeline is unchanging, the dynamic transporting model of natural gas pipeline is [45]

(3)\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }\left(\rho v\right)}{\mathrm{\partial }x}}=0,\\ {\displaystyle \frac{\mathrm{\partial }\left(\rho v\right)}{\mathrm{\partial }t}}+{\displaystyle \frac{\mathrm{\partial }\left(\rho v^2+p\right)}{\mathrm{\partial }x}}+{\displaystyle \frac{\lambda \rho v\left|v\right|}{2D}}=0,\end{array}

where p is gas pressure, ρ is gas density, v is gas velocity, \lambda is the Darcy friction factor, and D is the diameter of gas pipeline. The mass flow is q=\rho vS, where S denotes the section area.

Along with the actual state equation [17]

(4)\{\begin{array}{l}p=\rho Z{R}_{\mathrm{g}}T,\\ R_{\mathrm{g}}={\displaystyle \frac{R}{M}},\end{array}

where Z is compressibility factor, R_g is specific gas constant for mixed gas, R is general gas constant, T is gas temperature, and M is the molecular weight of the mixed gas.

Considering the hydrogen-blending situation, the change of \lambda needs to be studied due to the gradually mixing process along the pipeline. When the components of hydrogen-blending natural gas with is dynamically changing, the Darcy friction coefficient is different along the time and space coordinates. In this paper, hydraulic smooth and mixed friction zone are regarded as the two dominated influencing patterns in the natural gas flow [38], where the Darcy friction coefficient \lambda is calculated by Reynolds number Re and roughness \epsilon, i.e., \lambda =f(Re,\epsilon ). With Re>5\times {10}^{-4} and \epsilon >4\times {10}^{-4}, \lambda is calculated by

(5)\lambda =0.067{\left(2\epsilon +\frac{158}{Re}\right)}^{0.2}.

2.1.2 Model of compressor

The gas compressor in a natural gas pipeline network system is the main scheduling equipment, driven by the power of gas combustion or electricity [46].

(1) Electricity-driven compressor

The efficiency {\eta }_{\mathrm{e}\mathrm{d}} of the electric-driven compressor is calculated by the fitting the rectifying coefficients A_{\mathrm{e}\mathrm{d}} and B_{\mathrm{e}\mathrm{d}} with the real-world operating data [47] as

(6){\eta }_{\mathrm{e}\mathrm{d}}=A_{\mathrm{e}\mathrm{d}}{n}_{\mathrm{e}\mathrm{d}}+B_{\mathrm{e}\mathrm{d}},

where n_{\mathrm{e}\mathrm{d}} is the non-dimensional factor representing the operating condition of the compressor, defined as

(7)n_{\mathrm{e}\mathrm{d}}=\frac{P_{\mathrm{e}\mathrm{d}}}{P_{\mathrm{d}\mathrm{p},\mathrm{e}\mathrm{d}}},

where P_{\mathrm{e}\mathrm{d}} is the power consumption of the electricity-driven compressor and P_{\mathrm{d}\mathrm{p},\mathrm{e}\mathrm{d}} is the rated power in the condition of designing point.

With the definition of efficiency {\eta }_{\mathrm{e}\mathrm{d}} and the polytropic process working principles of the compressor, Eqs. (8) and (9) can be obtained.

(8)P_{\mathrm{e}\mathrm{d}}=\frac{W_{\mathrm{p}\mathrm{o}\mathrm{l}}}{{\eta }_{\mathrm{e}\mathrm{d}}}{q}_{\mathrm{v}},

(9)W_{\mathrm{p}\mathrm{o}\mathrm{l}}=\frac{m}{m-1}{p}_{\mathrm{i}\mathrm{n}}{v^{\prime }}_{\mathrm{i}\mathrm{n}}\left[{\left(\frac{p_{\mathrm{o}\mathrm{u}\mathrm{t}}}{p_{\mathrm{i}\mathrm{n}}}\right)}^{\frac{m-1}{m}}-1\right],

where W_{\mathrm{p}\mathrm{o}\mathrm{l}} represents the effective power of the polytropic compressing process for per unit of gas volume, q_{\mathrm{v}} is the volume flow of natural gas through the compressor and its adjacent pipeline inlet and outlet nodes, {v^{\prime }}_{\mathrm{i}\mathrm{n}} is the specific volume ratio of inlet natural gas, m is the polytropic coefficient of natural gas, and p_{\mathrm{o}\mathrm{u}\mathrm{t}} and p_{\mathrm{i}\mathrm{n}} represent the compressor outlet and inlet pressure, respectively [48].

(2) Gas combustion-driven compressor

Similar to the electric-driven compressor, the efficiency {\eta }_{\mathrm{g}\mathrm{d}} of the gas combustion-driven compressor is obtained by fitting the rectifying coefficients A_{\mathrm{g}\mathrm{d}} and B_{\mathrm{g}\mathrm{d}} with the real-world operating data, as

(10){\eta }_{\mathrm{g}\mathrm{d}}=A_{\mathrm{g}\mathrm{d}}{n}_{\mathrm{g}\mathrm{d}}+B_{\mathrm{g}\mathrm{d}},

where n_{\mathrm{g}\mathrm{d}} is the non-dimensional factor representing the operating condition of the gas combustion-driven compressor, defined as

(11)n_{\mathrm{g}\mathrm{d}}=\frac{P_{\mathrm{g}\mathrm{d}}}{P_{\mathrm{d}\mathrm{p},\mathrm{g}\mathrm{d}}},

where P_{\mathrm{g}\mathrm{d}} is the power consumption of the gas combustion-driven compressor and P_{\mathrm{d}\mathrm{p},\mathrm{g}\mathrm{d}} is the rated power in the condition of designing point.

The power of the gas combustion-driven compressor satisfies Eqs. (12) and (13) [42].

(12)P_{\mathrm{g}\mathrm{d}}=\frac{W_{\mathrm{p}\mathrm{o}\mathrm{l}}}{{\eta }_{\mathrm{g}\mathrm{d}}}{q}_{\mathrm{v}},

(13)Q_{\mathrm{g}\mathrm{d}}=\frac{P_{\mathrm{g}\mathrm{d}}}{{\eta }_{\mathrm{g}}{H}_{\mathrm{m}\mathrm{i}\mathrm{x}}},

where Q_{\mathrm{g}\mathrm{d}} is the volume flow of the natural gas cost of the gas combustion-driven compressor, {\eta }_{\mathrm{g}} is the efficiency of gas combustion, and H_{\mathrm{m}\mathrm{i}\mathrm{x}} is the calorific value of the natural gas–hydrogen mixed gas.

2.2 Modeling of power grid

2.2.1 Active current power flow model of power system

Power flow refers to the process in which the current or power flows from the power supply to the user load through various components of the system under the excitation of the power supply potential during the operation of the power system. Power flow calculation is to determine the operating state of the system under given operating conditions, such as voltage (amplitude and phase angle) of each bus, power distribution and power loss in the network. The steady-state model of the power system in IES adopts the active current (AC) power flow model, and the state variables of the system are mainly node voltage and phase angle. The balance equation of the active and reactive power injected by node is based on Ref. [49], as shown in Eq. (14).

(14)\{\begin{array}{c}{P}_i=U_i\sum _{j\in i}{U}_j\left(G_{ij}\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_i-{\theta }_j\right)+B_{ij}\mathrm{s}\mathrm{i}\mathrm{n}\left({\theta }_i-{\theta }_j\right)\right),\\ V_i=U_i\sum _{j\in i}{U}_j\left(G_{ij}\mathrm{s}\mathrm{i}\mathrm{n}\left({\theta }_i-{\theta }_j\right)-B_{ij}\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_i-{\theta }_j\right)\right),\end{array}

where P_i and V_i are the active and reactive power injected into Node i; U_i and U_j are the voltage amplitudes of Nodes i and j; G_{ij} and B_{ij} are the real and imaginary parts of the elements in Row i and Column j of the node admission matrix, respectively; and {\theta }_i and {\theta }_j are the phase angles of Nodes i and j, respectively.

Kirchhoff’s law in the circuit can be expressed as

(15)\left[\begin{array}{c}\mathit{A}\\ \mathit{B}\mathit{Z}\end{array}\right]\times \dot{\mathit{J}}=\left[\begin{array}{c}{\dot{\mathit{J}}}_{\mathbf{G}}-{\dot{\mathit{J}}}_{\mathbf{L}}\\ 0\end{array}\right],

where \mathit{A}\in {\mathit{R}}^{n\times b} is the correlation matrix of the branch, \mathit{B}\in {\mathit{R}}^{n\times b} is the basic loop matrix, \mathit{Z}\in {\mathit{R}}^{n\times b} is the impedance matrix of the branch, \dot{\mathit{J}}\in {\mathit{R}}^{\mathit{b}\times 1} is the current vector of the branch, {\dot{\mathit{J}}}_{\mathbf{G}}\in {\mathit{R}}^{n\times b} is current source vector and {\dot{\mathit{J}}}_{\mathbf{L}}\in {\mathit{R}}^{n\times b} is current load vector, n is the number of nodes, b is the number of branches, and k is the number of independent circuits.

2.2.2 Model of generator

1) Coal-fired generator

Coal-fired units use the heat released by coal to transfer heated water to steam for driving a turbine, which then drives a generator to generate electricity. Coal quantity M_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}} and power generation P_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}} satisfy Eq. (16).

(16)P_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}}=M_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}}{H}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}}{\eta }_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}},

where H_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}} is the calorific value of coal combustion, and {\eta }_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}} is the power generation efficiency of the coal-fired generator.

2) Gas turbine generator

With gas turbine exhaust high temperature and pressure gas into steam, the steam turbine is driven to generate electricity. As one of the coupled devices of the power grid and gas network, the calculation formula of power generation P_{\mathrm{g}\mathrm{a}\mathrm{s}} is

(17)P_{\mathrm{g}\mathrm{a}\mathrm{s}}=Q_{\mathrm{g}\mathrm{a}\mathrm{s}}{H}_{\mathrm{g}\mathrm{a}\mathrm{s}}{\eta }_{\mathrm{g}\mathrm{a}\mathrm{s}},

where Q_{\mathrm{g}\mathrm{a}\mathrm{s}} is gas consumption flow, H_{\mathrm{g}\mathrm{a}\mathrm{s}} is the calorific value of gas, and {\eta }_{\mathrm{g}\mathrm{a}\mathrm{s}} is the power generation efficiency of the gas unit.

(3) Renewable energy plants and water electrolysis hydrogen production station

Due to the diversity of renewable energy, solar, wind and water energy are considered in this paper, and the power generation of renewable energy plants is simulated according to local solar irradiance, wind speed, water flow and other meteorological parameters. The power generated by renewable energy sources is unstable, marked by P_{\mathrm{r}\mathrm{e}}. A part of P_{\mathrm{r}\mathrm{e}} is utilized and consumed by the power grid, marked with P_{\mathrm{r}\mathrm{e},\mathrm{u}}, while the other part is surplus and further used to produce hydrogen by electrolysis of water, marked with P_{\mathrm{r}\mathrm{e},\mathrm{s}}.

(18)P_{\mathrm{r}\mathrm{e}}=P_{\mathrm{r}\mathrm{e},\mathrm{u}}+P_{\mathrm{r}\mathrm{e},\mathrm{s}}.

The mathematical description of the hydrogen production process is

(19)Q_{{\mathrm{H}}_2}=\frac{P_{\mathrm{r}\mathrm{e},\mathrm{s}}{\eta }_{{\mathrm{H}}_2}}{H_{{\mathrm{H}}_2}},

where Q_{{\mathrm{H}}_2} is the standard volume flow rate of the hydrogen produced, P_{\mathrm{r}\mathrm{e}} is the excess power of renewable energy, and {\eta }_{{\mathrm{H}}_2} denotes the efficiency of the hydrogen production process.

3 Methodology

The two key issues of dynamic simulation and characteristics analysis for GEH-IES with hydrogen-blending transportation are to tackle the interaction of coupled variables between the power grid and gas networks, and to solve the transient flowing process through the hydrogen-blending natural gas pipeline. In this section, a novel dynamic simulation method for GEH-IES with hydrogen-blending is proposed to overcome the above problems. An iterative-interactive dynamic simulation method is presented to fulfill the overall solving framework of GEH-IES. Furthermore, based on a cell-segment method, an inner solving scheme for dynamic mixing process through the natural gas pipeline is proposed and integrated into the overall solving framework. The proposed method first enables the dynamic simulation of GEH-IES along with the distributing simulation of the hydrogen-blending gas networks.

3.1 Iterative-interactive dynamic simulation method for GEH-IES with hydrogen-blending transportation

A power grid is composed of distribution networks, user nodes, and source nodes. The user nodes include the electricity-driven compressor unit in the gas network in addition to the general power grid user nodes. The source nodes include gas-fired power plants, coal-fired power plants, and renewable power plants, of which the gas-fired power plants is also the user node of the gas network. Due to the instability of the renewable energy, the surplus electricity is used for hydrogen production, which is further mixed into the hydrogen mixing node of the gas network, affecting the flow dynamic characteristics of the whole gas pipeline and the performance of the compressor units. In this paper, natural gas network mainly includes gas transmission pipeline, gas/electricity-driven compressor units, and user nodes. As the user nodes of the power grid, the electricity-driven compressors affect all node states of the power grid. The detailed coupled relationship of the two systems in a GEH-IES is demonstrated in Fig.1.

Since the operation of hydrogen-blending complicates the dynamics of the GEH-IES, the time-spatial distribution of the natural gas pipeline networks and the interactive dynamics of hydrogen-blending nodes coupled between the two systems is urgently to be solved. In this section, an iterative-interactive dynamic simulation method for such a GEH-IES is proposed to tackle the problems, as shown in Fig.1. The solving scheme of the proposed method is decomposed into iterative calculating steps, of which the state parameters of the gas network and the power grid are interacting in each time step. One simulating step is further decomposed into three main tasks, updating the parameters of power grid, calculating parameters of hydrogen-blending node, and updating the parameters of gas network, as shown in Fig.2.

Considering the need of the distributing simulation for the natural gas pipeline, the overall time step is set as \mathrm{\Delta }t=60s. The other reason for such fine-grained time resolution is to satisfy the solving precision when the parameters of coupled nodes is one-step-time delayed, which is to make a closure of the iterative solving framework.

The presented solving algorithm for the GEH-IES with hydrogen-blending transportation is conducted as the follows. For the purpose of expressing the solving scheme clearly, several parameter vectors at time step t are first defined. The parameter vectors of the power grid include voltage amplitude {\mathit{U}}^t, generator active power {\mathit{P}}_{\mathrm{g}}, and load demand power {\mathit{P}}_{\mathrm{d}}^t:

{\mathit{U}}^t={(U_0^t,U_1^t,\dots ,U_i^t,\dots )}^{\mathrm{T}},

{\mathit{P}}_{\mathrm{g}}^t={(P_{0,\mathrm{g}}^t,P_{1,\mathrm{g}}^t,\dots ,P_{i,\mathrm{g}}^t,\dots )}^{\mathrm{T}}={(P_{\mathrm{b}\mathrm{l}\mathrm{c}}^t,{{\mathit{P}}_{\mathrm{g}\mathrm{a}\mathrm{s}}^t}^{\mathrm{T}},{{\mathit{P}}_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{l}}^t}^{\mathrm{T}},{{\mathit{P}}_{\mathrm{r}\mathrm{e},\mathrm{s}}^t}^{\mathrm{T}},{{\mathit{P}}_{\mathrm{r}\mathrm{e},\mathrm{u}}^t}^{\mathrm{T}})}^{\mathrm{T}},

(20){\mathit{P}}_{\mathrm{d}}^t={(P_{0,\mathrm{d}}^t,P_{1,\mathrm{d}}^t,\dots ,P_{i,\mathrm{d}}^t,\dots )}^{\mathrm{T}}={({{\mathit{P}}_{\mathrm{e}\mathrm{d}}^t}^{\mathrm{T}},{{\mathit{P}}_{\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{r}}^t}^{\mathrm{T}})}^{\mathrm{T}},

where P_{\mathrm{b}\mathrm{l}\mathrm{c}}^t is the active power of balance-generator node, and {\mathit{P}}_{\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{r}}^t is the vector of the users’ load demand power except the electricity-driven compressors.

The parameter vectors of the gas network include mass flow q^t, gas pressure p^t, hydrogen mixing ratio of H₂ r^t, and the gas loadQ^t (supply and demand):

{\mathit{q}}^t={(q_0^t,q_1^t,\dots ,q_i^t,\dots )}^{\mathrm{T}},

{\mathit{p}}^t={(p_0^t,p_1^t,\dots ,p_i^t,\dots )}^{\mathrm{T}},

{\mathit{r}}^t={(r_0^t,r_1^t,\dots ,r_i^t,\dots )}^{\mathrm{T}},

(21){\mathit{Q}}^t={(Q_0^t,Q_1^t,\dots ,Q_i^t,\dots )}^{\mathrm{T}}={({\mathit{Q}}_{\mathrm{g}\mathrm{a}\mathrm{s}}^t,{\mathit{Q}}_{\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{r}}^t,{\mathit{Q}}_{{\mathrm{H}}_2}^t)}^{\mathrm{T}},

where {\mathit{Q}}_{\mathrm{g}\mathrm{a}\mathrm{s}}^t is the demand vector of the gas consumption of gas-fired power plants, {\mathit{Q}}_{\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{r}}^t is the demand vector of users in the gas network except gas-fired power plants, and {\mathit{Q}}_{{\mathrm{H}}_2}^t is the supply vector of H₂. On the hydrogen injection node, the hydrogen mixing ratio r is calculated by

(22)r=\frac{Q_{{\mathrm{H}}_2}{\rho }_{{\mathrm{H}}_2}}{Q_{{\mathrm{H}}_2}{\rho }_{{\mathrm{H}}_2}+q},

where {\rho }_{{\mathrm{H}}_2} is the density of hydrogen in the standard state. Given the dynamic boundary condition of the two systems, {\mathit{p}}_{\mathrm{b}}\left(t\right), {\mathit{Q}}_{\mathrm{b}}\left(t\right), {\mathit{Q}}_{\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{r}}^{}\left(t\right), {\mathit{P}}_{\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{r}}^{}\left(t\right), and {\mathit{P}}_{\mathrm{r}\mathrm{e}}^{}\left(t\right), the proposed iterative-interactive simulation method for the GEH-IES is conducted as follows. Note that {\mathit{p}}_{\mathrm{b}}\left(t\right) and {\mathit{Q}}_{\mathrm{b}}\left(t\right) are the functions of inlet pressure and outlet mass flow on the boundary nodes in the gas network, respectively.

(1) Initialize the state parameters of the power grid and the gas network based on the results of static simulation of the GEH-IES at t=0s.

The vectors of state variables of power grid, {\mathit{U}}^t, {\mathit{\theta }}^t, and {\mathit{P}}_{\mathrm{g}}^t are determined via the optimal flow solution acquired by minimizing the reactive power:

\mathrm{m}\mathrm{i}\mathrm{n}F={\sum }_{i=1}^{N_{\mathrm{g}}}{P}_{i,\mathrm{g}}^t-{\sum }_{i=1}^{N_{\mathrm{d}}}{P}_{i,\mathrm{d}}^t,

(23)

where N_{\mathrm{g}} and N_{\mathrm{d}} are the number of generator and demand nodes in the power grid, respectively. Since {\mathit{P}}_{\mathrm{g}\mathrm{a}\mathrm{s}}^t in {\mathit{P}}_{\mathrm{g}}^t and {\mathit{P}}_{\mathrm{e}\mathrm{d}}^t in {\mathit{P}}_{\mathrm{d}}^t is of coupled nodes with the gas network, the extra constrains of the optimization problem in Eq. (23) are considered as the coupled relation of the gas-fired power plant power {\mathit{P}}_{\mathrm{g}\mathrm{a}\mathrm{s}}^t and gas demand {\mathit{Q}}_{\mathrm{g}\mathrm{a}\mathrm{s}}^t (see Eq. (17)), the coupled relation of the power of the electricity-driven compressor {\mathit{P}}_{\mathrm{e}\mathrm{d}}^t and the parameters of the gas networks (see Eqs. (6)–(9)).

The static solution to the gas network is determined by the model of the steady-state gas network and compressors (see Eqs. (1) and (2) and Eqs. (6)–(13)).

(2) Conduct one iterative step, \mathit{t}\leftarrow \mathit{t}+\mathrm{\Delta }\mathit{t}. The hydrogen-blending operation is executed and calculated.

Instead of producing and injecting hydrogen into the gas network presently, hydrogen-blending operation is one-time-step delayed to make a closure for the iterative scheme, which is much closer to the real situation. According to the solution to the generator power vector {\mathit{P}}_{\mathrm{g}}^{t-\mathrm{\Delta }t}, i.e., {\mathit{P}}_{\mathrm{r}\mathrm{e},\mathrm{u}}^{t-\mathrm{\Delta }t}, the volume flow {\mathit{Q}}_{{\mathrm{H}}_2}^{t-\mathrm{\Delta }t} of the injected H₂ is calculated via Eq. (19).

Likewise, the coupled variables {\mathit{Q}}_{\mathrm{g}\mathrm{a}\mathrm{s}}^t and {\mathit{P}}_{\mathrm{e}\mathrm{d}}^t is one-time-step delayed for the gas network simulation, too.

(3) Update the parameters of the gas network at time t.

According to the boundary condition functions {\mathit{p}}_{\mathrm{b}}\left(t\right), {\mathit{Q}}_{\mathrm{b}}\left(t\right), and {\mathit{Q}}_{\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{r}}^{}\left(t\right), along with the solved injection volume flow {\mathit{Q}}_{{\mathrm{H}}_2}^{t-\mathrm{\Delta }t} and load demand {\mathit{Q}}_{\mathrm{g}\mathrm{a}\mathrm{s}}^{t-\mathrm{\Delta }t}, the inner distribution simulation of the hydrogen-blending natural gas flow is conducted, which will be described in Section 3.2. The parameters such as the mass flow {\mathit{q}}^t, gas pressure {\mathit{p}}^t, and hydrogen mixing ratio of H₂ {\mathit{r}}^t are solved.

(4) Update the parameters of the power grid at time t.

Calculate the {\mathit{P}}_{\mathrm{e}\mathrm{d}}^t via Eqs. (6)–(9). According to the boundary condition functions {\mathit{P}}_{\mathrm{u}\mathrm{s}\mathrm{e}\mathrm{r}}^{}\left(t\right) and {\mathit{P}}_{\mathrm{r}\mathrm{e}}^{}\left(t\right), acquire the parameters of the power grid by solving Eq. (22) and update {\mathit{P}}_{\mathrm{g}}^t. Then back to (2).

3.2 Inner solving scheme for dynamic mixing process in natural gas pipeline

Since the operation of hydrogen blending is unstable due to the fluctuation of renewable energy, the mixing process of hydrogen with natural gas has a great impact on the property of such mixed gas, and further influences transporting dynamics of the pipeline flow. In this section, a cell-segmented method is proposed to capture the dynamics of hydrogen-blending transporting, fulfilling the dynamic distributing simulation for the gas pipeline. The fourth–fifth order Runge-Kutta algorithm [50] with adaptive step is utilized as the basic ODEs solver in this paper.

According to q=\rho vS and gas state Eq. (4), Eq. (3) turns to

(24)\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }t}}=-{\displaystyle \frac{Z{R}_{\mathrm{g}}T}{S}}{\displaystyle \frac{\mathrm{\partial }q}{\mathrm{\partial }x}},\\ {\displaystyle \frac{\mathrm{\partial }q}{\mathrm{\partial }t}}=-{\displaystyle \frac{Z{R}_{\mathrm{g}}T}{S}}\left({\displaystyle \frac{2q}{p}}{\displaystyle \frac{\mathrm{\partial }q}{\mathrm{\partial }x}}-{\displaystyle \frac{q^2}{p^2}}{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}}\right)-S{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}}-{\displaystyle \frac{\lambda Z{R}_{\mathrm{g}}T}{2SD}}{\displaystyle \frac{q^2}{p}}.\end{array}

To simulate the mixing process of the hydrogen within the natural gas, the gas pipeline is divided into multiple cells with given volume. In each cell, physical parameters like hydrogen mixing ratio and gas pressure are assumed to be independent of spatial distribution, meaning that in each cell only the ODEs problem is involved, and the mass flow rate is defined in the inlet and outlet interfaces of each cell. With the format of difference quotient, \mathrm{\Delta }q/\mathrm{\Delta }x and \mathrm{\Delta }p/\mathrm{\Delta }x, the dominating ODEs of the parameters in a cell can be written as

(25)\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}p}{\mathrm{\partial }t}}=-\alpha \mathrm{\Delta }q,\\ {\displaystyle \frac{\mathrm{d}q}{\mathrm{\partial }t}}=-\alpha {\displaystyle \frac{2q}{p}}\mathrm{\Delta }q+\left(\alpha {\displaystyle \frac{q^2}{p^2}}-{\displaystyle \frac{S}{\mathrm{\Delta }x}}\right)\mathrm{\Delta }p-\beta {\displaystyle \frac{q^2}{p}},\end{array}

(26)\alpha =\frac{Z{R}_{\mathrm{g}}T}{S\mathrm{\Delta }x},\beta =\frac{\lambda Z{R}_{\mathrm{g}}T}{2SD},

where α and β are determined by gas physical properties (compression factor Z, the gas constant R_g, and the hydraulic friction coefficient λ). Flow variables including pressure p of the cell and mass flow of the cell bound (or nodes) q are changed by tracking the gas components in real time.

With the gas flows, the parameters of adjacent cells are interacting with each other explicitly. To fulfill the cell-segmented solving scheme, the storage \sigma of the cell is the mass of gas, as shown in Fig.3, where {\sigma }_i^t represents storage of the ith cell at time t, which meets the conservation of mass as follows. Equation (26) provides a way to solve the ODEs system of cells time, separately.

(27){\sigma }_i^t={\sigma }_i^{t-\mathrm{\Delta }t}+\mathrm{\Delta }t{q}_i^t-\mathrm{\Delta }t{q}_{i+1}^t,

where q_i^t is the inlet mass flow of the ith cell at time t, corresponding to the nodes in the gas networks. The average density {\rho }_i^t of the cell can be calculated, as expressed in Eq. (28).

(28){\rho }_i^t=\frac{{\sigma }_i^t}{S\mathrm{\Delta }x}.

Supposing the components of natural gas is unchanged and the Z and R_{\mathrm{g}} can be calculated based on hydrogen ratio r, the mixing ratio r_i^t of the ith cell can be obtained by Eq. (29) [51].

(29)p_i^t={\rho }_i^{t}Z(p_i^t,T)R_{\mathrm{g}}\left(r_i^t\right)T.

According to the formulation proposed by American Natural Gas Association [52], the compression factor Z is calculated by

(30)Z(p_i^t,T)=1+0.257\frac{p_i^t}{p_{\mathrm{c}}}-0.533\frac{p_i^t}{p_{\mathrm{c}}}\frac{T_{\mathrm{c}}}{T},

where p_{\mathrm{c}} and T_{\mathrm{c}} are quasi-critical pressure and quasi-critical temperature evaluated according to component ratio, and R_{\mathrm{g}} is determined by the molecular weight of gas M,

(31)M=\frac{{M_{{\mathrm{H}}_2}M}_{\mathrm{n}\mathrm{g}}}{M_{\mathrm{n}\mathrm{g}}{r}_i^t+M_{{\mathrm{H}}_2}(1-r_i^t)}.

Thus,

(32)R_{\mathrm{g}}\left(r_i^t\right)=\frac{R(M_{\mathrm{n}\mathrm{g}}{r}_i^t+M_{{\mathrm{H}}_2}\left(1-r_i^t\right))}{{M_{{\mathrm{H}}_2}M}_{\mathrm{n}\mathrm{g}}},

where M_{{\mathrm{H}}_2} and M_{\mathrm{n}\mathrm{g}} are the molar mass of hydrogen and natural gas. The viscosity µ of mixed gas is calculated by Wilke’s semi-empirical relation, as Eq. (33) [53].

(33)\mu =\frac{{\mu }_{{\mathrm{H}}_2}{r}_i^t\sqrt{M_{{\mathrm{H}}_2}}+{\mu }_{\mathrm{n}\mathrm{g}}\left(1-r_i^t\right)\sqrt{M_{\mathrm{n}\mathrm{g}}}}{r_i^t\sqrt{M_{{\mathrm{H}}_2}}+\left(1-r_i^t\right)\sqrt{M_{\mathrm{n}\mathrm{g}}}},

where {\mu }_{{\mathrm{H}}_2} and {\mu }_{\mathrm{n}\mathrm{g}} are the viscosity of hydrogen and natural gas. Then Re is calculated by [38]

(34)Re=\frac{{\rho }_i^t{v}_i^{t}D}{\mu }=\frac{q_i^{t}D}{S\mu }.

Thus,

(35)\lambda \left(r_i^t\right)=0.067{\left(2\epsilon +\frac{158}{Re}\right)}^{0.2}.

Rewrite ODEs for the ith cell,

(36)\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{p}_i}{\mathrm{d}t}}=-{\alpha }_i\left(q_{i+1}-q_i\right),\\ {\displaystyle \frac{\mathrm{d}{q}_i}{\mathrm{d}t}}=-{\alpha }_i{\displaystyle \frac{2q_i}{p_i}}\left(q_{i+1}-q_i\right)+\left({\alpha }_i{\displaystyle \frac{{q_i}^2}{{p_i}^2}}-{\displaystyle \frac{S}{\mathrm{\Delta }x}}\right)\left(p_{i+1}-p_i\right)-{\beta }_i{\displaystyle \frac{{q_i}^2}{p_i}},\end{array}

(37){\alpha }_i=\frac{Z(p_i^t,T)R_{\mathrm{g}}\left(r_i^t\right)T}{S\mathrm{\Delta }x}{\beta }_i=\frac{\lambda \left(r_i^t\right)Z(p_i^t,T)R_{\mathrm{g}}\left(r_i^t\right)T}{2SD}.

Note that parameters α_i and β_i of the ith cell remain unchanged during \mathrm{\Delta }t. By solving the DAEs constructed by Eqs. (27)–(30), Eq. (32) and Eqs. (33)–(37) with i=0,1,\dots ,n, where number of the cells along the gas pipeline n, the parameters of gas network at time t, mass flow {\mathit{q}}^t, gas pressure p^t and concentration of H₂{\mathit{r}}^t are obtained. Tab.1 tabulates the overall flow of the proposed simulation framework.

4 Case study

4.1 Validation of dynamic hydrogen-mixing process solving method

The proposed inner solving scheme for dynamic mixing process in natural gas pipeline contributes significantly to the precision of flow parameters and hydrogen mixing ratio in the whole GEH-IES simulation. To validate the method, the model for a real-world hydrogen-blending natural gas pipeline is built, which is a coal-to-gas pipeline from Xinjiang Uygur Autonomous Region in China [35]. The hydrogen is produced by coal and then injected to the inlet node as shown in Fig.4. The structure parameters of the pipeline are listed in Tab.2.

The monitoring parameters is gas flow rate q_{\mathrm{i}\mathrm{n}}, gas pressure p_{\mathrm{i}\mathrm{n}}, hydrogen mole fraction in the inlet node, and the gas pressure p_{\mathrm{o}\mathrm{u}\mathrm{t}} in the outlet node. Here, the gas flow rate q_{\mathrm{i}\mathrm{n}} and hydrogen mole fraction on the inlet node and the gas pressure p_{\mathrm{o}\mathrm{u}\mathrm{t}} on the outlet node are used as the boundary conditions, as shown in Fig.5. The true value and simulated value of gas pressure p_{\mathrm{i}\mathrm{n}} on the inlet node is compared to validate the dynamic hydrogen-mixing process solving method in Section 3.2.

As shown in Fig.6, the calculated gas pressure p_{\mathrm{i}\mathrm{n}} precisely follows its true value, with an overall error less than 0.02%, proving the computational accuracy of the proposed inner solving scheme for dynamic mixing process in natural gas pipeline. Furthermore, the distribution of hydrogen mole fraction along the simulating time interval is shown as Fig.7 while the change of hydrogen mole fraction at different locations along the simulating time interval is shown in Fig.8. With the increase of time, the hydrogen concentration gradually approaches to stability.

4.2 Cases of real-world GEH-IES under two dynamic operating conditions

In this paper, two typical operating conditions, the daily fluctuating renewable energy and the sudden shutdown of electricity–power generator, are mainly simulated and analyzed based on the proposed method to understand the dynamics of the GEH-IES with hydrogen-blending transportation. A GEH-IES is built referring to the real-world system. The real natural gas network in a region of central China is selected as the gas network, and the coupled power grid is simplified as an IEEE 39-node grid [54]. The whole system structure is shown in Fig.9. The gas network can be decomposed into an upper and a lower road, from Nodes 1 to 3 and from Nodes 5 to 14, respectively. The parameters of the inlet and outlet nodes of the two road and the compressor performance are mainly analyzed in this section.

The power grid system consists of 10 power generation nodes. The natural gas network includes four compressor units, where Nodes 2, 5, and 13 are electricity-drive compressors while Node 9 is a gas combustion-driven compressor. Nodes 1, 4, 7, 8, 10, 11, and 12 are load nodes. The hydrogen produced by surplus renewable energy is mixed into Nodes 1, 7, and 8. For more details on the system structure, please refer to Tab.2 to Tab.4.

The components of natural gas are listed in Tab.6 [15]. The two main components are significantly different.

4.2.1 Daily fluctuation of renewable energy

The system is simulated under the condition of daily fluctuation of renewable energy power generation while the other boundary conditions are set as constant, and the daily fluctuation curve of renewable energy power generation is shown in Fig.10. Due to the lack of local photovoltaic resources, the power of photovoltaic power stations is less. Because high hydrogen concentration may lead to the occurrence of hydrogen embrittlement, making the performance of the pipeline steel degenerate, which may lead to risks such as leaks. Taking into account the safety requirements of the pipeline, the hydrogen mixing ratio after hydrogen-blending should not be higher than 10%. In addition, to avoid the effects of the initial condition, the GEH-IES is simulated for about 20 days. The simulation data from day 10 were used for analysis.

Fig.11 shows the outlet pressure and inlet mass flow response of the lower and upper road when the renewable energy fluctuates within the day. Since the amount of hydrogen injection to the upper road is relatively small, the effect of hydrogen mixture of outlet pressure is minimal on the upper road, as shown in Fig.11 (a) and Fig.11(c).

Moreover, Fig.11(c) and Fig.11(d) indicate that the fluctuation of the upper and lower road flow in a day basically depends on the fluctuation of the surplus renewable power used to produce hydrogen. Since the upper road is shorter, the inlet flow shows a stepped and fast response. The lower road is longer, resulting in a smoother inlet flow curve with an obvious time delay. Thus, it can be concluded that the variation of the outlet pressure is basically opposite to that of the inlet flow when the power of renewable generators changes, and both of the responses are time-delayed to different extent in respect to the distance of mixing flow. It is worth noting that the flow mass is more sensitive than the pressure.

Fig.12 shows the change of consumption of four compressor units when the power of renewable generators changes. Fig.12(a) and Fig.12(c) show that the compressor energy consumption of Node 2 basically changes with the upper road flow. A comparison of Fig.12(a) and Fig.12(b) with Fig.12(c) and Fig.12(d) suggests that the closer the compressor unit to the hydrogen blending node is, where the hydrogen is more likely not to be evenly mixed, the larger the power fluctuation is, proving that the hydrogen mixing state has a great impact on the performance of the compressor. The power of the electricity-driven compressor is mainly influenced by flow parameters, while the gas combustion-driven compressor is depending on both flow parameters and gas components. When the gas combustion-driven compressor is closer to the hydrogen blending node, stepped response can be observed.

Fig.13 shows the response of power and voltage of nodes considering renewable energy daily fluctuation. Due to the fact that there are a large number of power grid nodes, only one curve is shown for each type of generation type. It can be seen from the results that the power of renewable energy generator fluctuates greatly, which is mainly caused by the volatility of renewable energy. The power of coal-fired generator also has certain volatility, because coal-fired generators need to undertake the task of peak load. Fig.13(b) shows the voltage amplitude of node relative to the rated voltage. As can be seen from Fig.13(b), the voltage fluctuates around the rated voltage, and the deviation is not more than 6%.

Fig.14 shows the change of hydrogen concentration of different roads at the output node when the surplus power of renewable generators fluctuates within the day. The hydrogen concentration of output mainly changes with the hydrogen production near the inlet node.

Fig.15 shows the response of outlet density and outlet velocity of the upper and lower roads considering renewable energy daily fluctuation. The results demonstrates that the density curve is very similar to the pressure curve, because the density is affected by the pressure and the hydrogen mixture ratio, but the hydrogen mixture ratio changes little, and the main density is affected by the pressure. In addition, Fig.15(b) shows that the velocity of the upper road is much lower than that of the lower road, mainly because the flow rate of the road is small and the pressure is large, and only a small velocity is needed to meet the transport demand.

4.2.2 Sudden shutdown of generator

Assuming that the coal-fired generator at Node 36 suddenly stops at t = {10}^6 s while other boundary conditions of the system do not change, the response the flow parameters of GEH-IES with hydrogen-blending is compared with that of the traditional GE IES without hydrogen-blending.

Fig.16 reflects the change of hydrogen concentration at different roads at the outlet node when the generator shuts down. After the shutdown of the generator, the renewable energy generator needs to provide more energy for the power grid. Therefore, the hydrogen concentration of the system is reduced from 7.49% to 0.

Fig.17 shows the change of pressure at the output and flow at the input of different roads when the generator shuts down with and without hydrogen mixing. Since there is no generator load on the upper road, as shown in Fig.17(a) and 17(c), the upper road pressure and flow remain unchanged. As shown in Fig.17(b), on the lower road, the pressure at the output with hydrogen blending response is lower than that without hydrogen blending at the early stage since more hydrogen blending leads to more pressure loss, while both of the two cases reach the same pressure finally as the hydrogen mixing amount drops to 0. However, the outlet pressure with hydrogen blending first decreases and then increases, while the outlet pressure is decreasing all the time when there is no hydrogen blending. This is caused by the joint influence of hydrogen production and the gas flow. The shutdown of the generator leads to less hydrogen injection and more gas consumption of the gas-fired generator. Less hydrogen injection results in a less pressure loss due to the less molar mass of hydrogen, while more gas consumption results in a larger source gas mass flow. According to Fig.17(b), the gas consumption load responds quickly, while the respond to the change of hydrogen injection is relatively delayed.

A comparison of the two curves in Fig.17(d) suggests that the hydrogen injection leads to a lower demand of gas and a lower mass flow. As the amount of hydrogen injection decreases, the gas flow increases and then drops to the same level as that of no-hydrogen-involved condition. Nevertheless, the corresponding variables in the GE-IES without hydrogen injection is only affected by an increase in mass flow change.

Fig.18 shows the response of outlet density and outlet velocity of the lower road when the generator suddenly shuts down. As can be seen from Fig.18, the change in density and velocity with hydrogen injection is much greater than that without hydrogen injection. This is mainly because when the generator suddenly shuts down, in the case with hydrogen injection, the hydrogen concentration decreases rapidly, and the density and velocity are simultaneously affected by the change of pressure and hydrogen concentration. In the case without hydrogen injection, the density and velocity are only affected by pressure.

Fig.19 shows the response of the energy consumption of four compressor units when the generator shuts down with and without hydrogen injection. Since there is no generator load on the upper road, the compressors at the upper Node 2 remain unchanged with or without hydrogen injection, as shown in Fig.19(a). Fig.19(b) shows that the power curve of the compressor at Node 5 is positively correlated with the inlet flow mass (Fig.17(d)), since Node 5 is near the input of the lower pipeline, and the energy consumption of the compressor is mainly determined by the inlet mass flow.

As shown in Fig.19(c), since the gas combustion-driven compressor is located behind the hydrogen injection node, it is mainly impacted by the hydrogen mixing ratio and mass flow rate. This is shown in Fig.19(c) and Fig.19(d), since hydrogen blending results in a lower density of the gas flow, the flow volume increases and thus more power from gas are needed. With the shutdown of the gas-fired generator, the hydrogen mixing ratio decreases. Thus, the power of compressors at Nodes 9 and 13 starts to drop to a certain value smaller than that of the no-hydrogen-involved condition and then recovers to the same level.

Fig.20 shows the response of power and voltage of nodes when one of the generators shuts down suddenly with and without hydrogen injection. As can be seen from Fig.20, when a generator shuts down, in order to meet the demand for electricity, the power of some other generators is increased instantaneously. In addition, the power of some generators does not change because these generators have reached their maximum power. In the case with hydrogen injection, the variation of the power of the generator remains unchanged, but the value is different from that without hydrogen injection. In addition, the shutdown of the generator also causes the voltage amplitude of the node to change, but hydrogen injection has little effect on the voltage change.

5 Conclusions

In this paper, a novel dynamic simulation method for the GEH-IES with hydrogen-blending transportation is proposed. For the purpose of capturing the time-spatial dynamics of such a GEH-IES, an iterative-interactive simulation scheme for GEH-IES and a cell-segmented based inner solving method for hydrogen–natural gas mixing flow are proposed and fulfilled. A case based on the real-world system is built and simulated. To further disclose the dynamics of the GEH-IES, two typical dynamic conditions, surplus renewable energy fluctuation and power generator shutdown, are calculated and analyzed. Evolution patterns under certain settings of key parameters including outlet pressure, inlet mass flow, hydrogen mixing ratio and compressor power/gas consumption are mainly discussed. The dynamic characteristics of the GEH-IES are to some extent disclosed and the following conclusions can be drawn.

1) The proposed iterative-interactive solving scheme with the cell-segmented gas flow method makes it possible for the calculation of hydrogen mixing ratio in fine-grained time and space resolution, providing a dynamic time-spatial simulation and analysis method for GEH-IES.

2) The impact of hydrogen blending on flow parameters is to some extent time-delayed in respect of mixing distance, and further influences the performance of the compressor due to mixing degree. The more even of mixture, the less instability of flow parameters and the compressor power.

3) With the change of renewable energy in the power grid, the hydrogen mixing ratio and gas pressure loss along the gas pipeline are negative correlated. Along with the change of the mass flow, the gas pressure and the compressor power always vary a great deal then recover with a limited degree.

However, in order to ensure the safe and efficient operation of the GEH-IES, the control of the concentration of hydrogen and the joint optimization of the GEH-IES system need to be further explored.

References

[1] Wang C, Lv C, Li P. . Modeling and optimal operation of community IESs: A case study from China. Applied Energy, 2018, 230: 1242–1254

[2] Qin C, Yan Q, He G. IESs planning with electricity, heat and gas using particle swarm optimization. Energy, 2019, 188: 116044

[3] Huang D, Zhou D, Jia X. . A mixed integer optimization method with double penalties for the complete consumption of renewable energy in distributed energy systems. Sustainable Energy Technologies and Assessments, 2022, 52: 102061

[4] Sinsel S R, Riemke R L, Hoffmann V H. Challenges and solution technologies for the integration of variable renewable energy sources—A review. Renewable Energy, 2020, 145: 2271–2285

[5] ReynoldsRSlagerW. Pipeline transportation of hydrogen. In: Nejat Veziroğlu T, ed. Hydrogen Energy. New York: Springer, 1975, 533–543

[6] Mohammad A K, Sumeray C, Richmond M. . Assessing the sustainability of liquid hydrogen for future hypersonic aerospace flight. Aerospace, 2022, 9(12): 801

[7] Meng B, Guo F, Hu L. . Wind abandonment analysis of multi-energy systems considering gas–electricity coupling. Electric Power Engineering Technology, 2019, 38(6): 2–8 (in Chinese)

[8] Haeseldonckx D, D’haeseleer W. The use of the natural gas pipeline infrastructure for hydrogen transport in a changing market structure. International Journal of Hydrogen Energy, 2007, 32(10–11): 1381–1386

[9] Hanley E S, Deane J, Gallachóir B Ó. The role of hydrogen in low carbon energy futures—A review of existing perspectives. Renewable & Sustainable Energy Reviews, 2018, 82: 3027–3045

[10] Smit R, Weeda M, De Groot A. Hydrogen infrastructure development in The Netherlands. International Journal of Hydrogen Energy, 2007, 32(10–11): 1387–1395

[11] Saedi I, Mhanna S, Mancarella P. Integrated electricity and gas system modelling with hydrogen injections and gas composition tracking. Applied Energy, 2021, 303: 117598

[12] Zhao J, Chen L, Wang Y. . A review of system modeling, assessment and operational optimization for IESs. Science China. Information Sciences, 2021, 64(9): 1–23

[13] Shen F, Ju P, Shahidehpour M. . Singular perturbation for the dynamic modeling of IESs. IEEE Transactions on Power Systems, 2020, 35(3): 1718–1728

[14] Ma S, Sun S, Wu H. . Decoupling optimization of IES based on energy quality character. Frontiers in Energy, 2018, 12(4): 540–549

[15] Tabkhi F, Azzaro-Pantel C, Pibouleau L. . A mathematical framework for modelling and evaluating natural gas pipeline networks under hydrogen injection. International Journal of Hydrogen Energy, 2008, 33(21): 6222–6231

[16] Agaie B, Khan I, Yacoob Z. . A novel technique of reduce order modelling without static correction for transient flow of non-isothermal hydrogen–natural gas mixture. Results in Physics, 2018, 10: 532–540

[17] Guandalini G, Colbertaldo P, Campanari S. Dynamic modeling of natural gas quality within transport pipelines in presence of hydrogen injections. Applied Energy, 2017, 185: 1712–1723

[18] Guandalini G, Colbertaldo P, Campanari S. Dynamic quality tracking of natural gas and hydrogen mixture in a portion of natural gas grid. Energy Procedia, 2015, 75: 1037–1043

[19] Cheli L, Guzzo G, Adolfo D. . Steady-state analysis of a natural gas distribution network with hydrogen injection to absorb excess renewable electricity. International Journal of Hydrogen Energy, 2021, 46(50): 25562–25577

[20] Elaoud S, Hadj-Taïeb E. Transient flow in pipelines of high-pressure hydrogen–natural gas mixtures. International Journal of Hydrogen Energy, 2008, 33(18): 4824–4832

[21] Liu B, Liu S, Guo S. . Economic study of a large-scale renewable hydrogen application utilizing surplus renewable energy and natural gas pipeline transportation in China. International Journal of Hydrogen Energy, 2020, 45(3): 1385–1398

[22] ITMPower. HyDeploy: UK gas grid injection of hydrogen in full operation. 2020-6-1, available at the website of ITM Power

[23] Suzuki T, Kawabata S, Tomita T. Present status of hydrogen transport systems utilizing existing natural gas supply infrastructures in Europe and the USA. The Institute of Energy Economics (Japan) Report, 2006

[24] TiekstraGKoopmanF. The NATURALHY project: First step in assessing the potential of the existing natural gas network for hydrogen delivery. In: International Gas Union Research Conference, Paris, France, 2008

[25] Colbertaldo P, Guandalini G, Campanari S. Modelling the integrated power and transport energy system: The role of power-to-gas and hydrogen in long-term scenarios for Italy. Energy, 2018, 154: 592–601

[26] Fang R. Life cycle cost assessment of wind power–hydrogen coupled IES. International Journal of Hydrogen Energy, 2019, 44(56): 29399–29408

[27] Pan G, Gu W, Qiu H. . Bi-level mixed-integer planning for electricity–hydrogen IES considering levelized cost of hydrogen. Applied Energy, 2020, 270: 115176

[28] He J, Wu Y, Wu M. . Two-stage configuration optimization of a novel standalone renewable IES coupled with hydrogen refueling. Energy Conversion and Management, 2022, 251: 114953

[29] Fang R. Multi-objective optimized operation of IES with hydrogen storage. International Journal of Hydrogen Energy, 2019, 44(56): 29409–29417

[30] Li N, Zhao X, Shi X. . IESs with CCHP and hydrogen supply: A new outlet for curtailed wind power. Applied Energy, 2021, 303: 117619

[31] Hajimiragha A, Canizares C, Fowler M, et al. Optimal energy flow of IESs with hydrogen economy considerations. In: 2007 iREP Symposium—Bulk Power System Dynamics and Control-VII. Revitalizing Operational Reliability. Charleston, SC, USA, 2007

[32] Liu J, Xu Z, Wu J. . Optimal planning of distributed hydrogen-based multi-energy systems. Applied Energy, 2021, 281: 116107

[33] Zhou D, Yan S, Huang D. . Modeling and simulation of the hydrogen blended gas−electricity IES and influence analysis of hydrogen blending modes. Energy, 2022, 239: 121629

[34] Liu J, Sun W, Harrison G P. The economic and environmental impact of power to hydrogen/power to methane facilities on hybrid power–natural gas energy systems. International Journal of Hydrogen Energy, 2020, 45(39): 20200–20209

[35] Chen X, Ma L, Zhou C. . Improved resistance to hydrogen environment embrittlement of warm-deformed 304 austenitic stainless steel in high-pressure hydrogen atmosphere. Corrosion Science, 2019, 148: 159–170

[36] Zhou D, Li T, Huang D. . The experiment study to assess the impact of hydrogen blended natural gas on the tensile properties and damage mechanism of X80 pipeline steel. International Journal of Hydrogen Energy, 2021, 46(10): 7402–7414

[37] Zhou C, Ye B, Song Y. . Effects of internal hydrogen and surface-absorbed hydrogen on the hydrogen embrittlement of X80 pipeline steel. International Journal of Hydrogen Energy, 2019, 44(40): 22547–22558

[38] Hardie D, Charles E, Lopez A. Hydrogen embrittlement of high strength pipeline steels. Corrosion Science, 2006, 48(12): 4378–4385

[39] WeiTZhouDYaoQ, . Simulation model uncertainty quantification and model calibration for natural gas compressor units. In: ASME International Mechanical Engineering Congress and Exposition, Salt Lake City, Utah, USA, 2019

[40] MaSZhouDZhangH, . Modeling and optimal operation of a network of energy hubs system with distributed energy resources. In: Turbo Expo: Power for Land, Sea, and Air, Charlotte, North Carolina, USA, 2017

[41] Wang C, Yan S, Shao R, et al. Coordinative optimization operation on the gas-electricity IESs. In: 2021 6th Asia Conference on Power and Electrical Engineering, Chongqing, China, 2021

[42] Yan S, Li T, Huang D, et al. Dynamic modeling and simulation of the large-scale regional integrated electricity and natural gas system. In: 2020 4th International Conference on Smart Grid and Smart Cities, Osaka, Japan, 2020

[43] YangDYanSZhouD, . Reinforcement learning methods on optimization problems of natural gas pipeline networks. In: 2020 4th International Conference on Smart Grid and Smart Cities, Osaka, Japan, 2020

[44] Liu C, Shahidehpour M, Wang J. Coordinated scheduling of electricity and natural gas infrastructures with a transient model for natural gas flow. Chaos, 2011, 21(2): 025102

[45] Kiuchi T. An implicit method for transient gas flows in pipe networks. International Journal of Heat and Fluid Flow, 1994, 15(5): 378–383

[46] Zhou D, Huang D, Hao J. . Vibration-based fault diagnosis of the natural gas compressor using adaptive stochastic resonance realized by generative adversarial networks. Engineering Failure Analysis, 2020, 116: 104759

[47] Zhou D, Wei T, Huang D. . Reliability assessment based on multisource information fusion method for high pressure natural gas compressors. Journal of Pressure Vessel Technology, 2021, 143(3): 031702

[48] Zhou D, Zhang H, Li Y G. . A dynamic reliability-centered maintenance analysis method for natural gas compressor station based on diagnostic and prognostic technology. Journal of Engineering for Gas Turbines and Power, 2016, 138(6): 061601

[49] Pan Z, Guo Q, Sun H. Interactions of district electricity and heating systems considering time-scale characteristics based on quasi-steady multi-energy flow. Applied Energy, 2016, 167: 230–243

[50] Dormand J R, Prince P J. A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics, 1980, 6(1): 19–26

[51] Liu K, Biegler L T, Zhang B. . Dynamic optimization of natural gas pipeline networks with demand and composition uncertainty. Chemical Engineering Science, 2020, 215: 115449

[52] MohringJHoffmannJHalfmannT, . Automated model reduction of complex gas pipeline networks. In: PSIG Annual Meeting, California, USA, 2004

[53] Bagajewicz M, Valtinson G. Computation of natural gas pipeline hydraulics. Industrial & Engineering Chemistry Research, 2014, 53(26): 10707–10720

[54] Zhou D, Ma S, Huang D. . An operating state estimation model for IESs based on distributed solution. Frontiers in Energy, 2020, 14(4): 801–816

RIGHTS & PERMISSIONS

Higher Education Press 2023