A Unique Electrical Model for the Steady-State Analysis of a Multi-Energy System

Abstract

In order to decarbonize the energy sector, the interdependencies between the power and natural gas systems are going to be much stronger in the next period. Thus, it is necessary to have a powerful simulation model that is able to efficiently and simultaneously solve all coupled energy carriers in a single simulation environment in only one simulation step. As an answer to the described computational challenges, a unique model for the steady-state analysis of a multi-energy system (MES) using the electrical analogy approach is developed. Detailed electrical equivalent models, developed using the network port theory and the load flow method formulation, of the most important natural gas network elements, as well as of the linking facilities between the power and natural gas systems, are given. The presented models were loaded up into a well-known software for the power system simulation—NEPLAN. In the case studies, the accuracy of the presented models is confirmed by the comparison of the simulation results with the results obtained by SIMONE—a well-known software for natural gas network simulations. Moreover, the applicability of the presented unique model is demonstrated by the MES security of a supply analysis.

1. Introduction

The main goal of actual energy policies is the decarbonization of the energy sector in the next medium to long term [1,2,3]. According to the joint Ten-Year Network Development Plan 2020 (TYNDP 2020) of ENTSOG and ENTSO-E [4], to comply with the 1.5 °C targets of the Paris Agreement [2], carbon neutrality in the power sector must be achieved by 2040, while in the overall energy sector by 2050. Consequently, the sector coupling between power and natural gas is recognized as the key contributor to achieving the decarbonization targets [5]. Accordingly, renewable energy sources (RES), primarily wind power plants (WPPs) and solar power plants (SPPs), will play an important role in the decarbonization of the power sector. Moreover, the decarbonization of the natural gas system will also have to play a significant role. For that purpose, a relatively new technology called power-to-gas (P2G) is recognized as one of the most important solutions [4,6]. The role of P2G in the decarbonization process is two-fold. On the one side, P2G technology will convert the surplus generation from intermittent RES into a so-called renewable gas, and help the power sector with variable RES integration [7], while, on the other side, P2G will help the gas sector to achieve the decarbonization targets [8,9]. Moreover, to balance the power system, when additional power injection is needed, e.g., due to a lower RES generation than predicted, the installation of new gas-fired power plants (GFPPs) will take part in the near future [4,10]. Moreover, the natural gas and power sectors are also already linked through the electrically driven gas compressors. In any case, the interdependency between these two energy infrastructures is going to be much stronger in the next period. Thus, the analyses of the coupled and strongly interdependent energy systems, as well as the operation and planning strategies, should be conducted on the integrated, so-called multi-energy system (MES) model [11].

MES has become the subject of research in numerous papers in the last 10–15 years [12,13,14]. From the beginning, interdependencies between the power and the natural gas systems have been mostly studied [15,16], which has continued to this day [17,18,19,20,21,22,23]. In [17,18], the influence of the gas network on solving the problems in the electric power network caused by wind power plants was analyzed. Models for the security of a supply analysis in integrated power and natural gas systems were developed in [19,20]. When we speak about the time resolution of the analysis of multi-energy systems, it can refer to the long-term planning of system development [21] or short-term management of the integrated system [22,23]. The main issues of all cited papers, as well as other papers dealing with the MES modeling, have been the complexity of modeling and slow computational performance. The former is a result of the applied modeling approach where every energy system is modeled separately and then linked through the consumption and constraints equations [24,25,26]. The latter is a direct consequence of the former, since, in that case, a high resource-consuming iterative inter-model calculation process must be used. For example, to conduct the security of a supply analysis in the power and natural gas networks coupled via GFPP [24], it is necessary to apply several calculation steps. In the first step, the level of GFPP engagement is determined by performing load flow calculations in the power network. After obtaining the GFPP generation, the associated gas consumption can be calculated in the second step. The GFPP gas consumption represents the gas load in the gas network, and the gas flow calculations in the gas network can be performed in the third step. However, if the overloading or outage of either network element occurs, that affects the GFPP operation, the described iterative inter-model calculations could be repeated several times before the calculation processes converge.

There are several software packages on the market that can be used to model and simulate more than one network (e.g., electrical, natural gas, water, and district heating networks). Unfortunately, this cannot be performed simultaneously in a single model, but in separate partial models. In [27], this issue is resolved with external managing by Python and the OPSIM module. However, both networks are still separated in their models, and the iterative inter-model process should also be applied.

The separate modeling of the coupled energy carrier infrastructures was also subject in [28,29,30]. In these papers, unlike the above-referenced papers, MES is analyzed simultaneously by using the Newton–Raphson approach, which has several significant limitations. The main constraints are a necessity to have the initial value of the iteration close enough to the solution, which can be an issue with the gas flow equations [29], and the system of equations has to be regular [31]. If initial gas pressures at the ends of the pipelines are close to each other, the Jacobi matrix becomes singular, and the results cannot be obtained by the Newton–Raphson technique. Moreover, every energy system is again separately described by its own set of equations that are only simultaneously solved if all Newton–Raphson pre-calculation conditions are satisfied. Consequently, every energy carrier still needs to have its slack node.

In this paper, the electrical analogy approach is applied in MES modeling that is consisted of the integrated power and natural gas systems. Hence, in addition to the elements of the power grid, the natural gas network elements, as well as the coupling facilities, such as P2G, GFPP, and electrically driven gas compressor stations, are in the unique MES model represented by the equivalent electrical elements. In the electrical analogy approach, the volumetric gas flow rate and gas pressure are analogous to the electric current and electric potential [32,33], respectively. Since, in this paper, the overall MES is modeled as a single power network, there is no more need for the above described high resource-consuming iterative inter-model calculations. Consequently, it is possible to simultaneously solve the coupled networks in only one step. Moreover, the applicability and simplicity of MES modeling are also achieved since the developed network element models can be easily used in the existing, robust, and well-known electrical software packages. The developed MES model is solved using the extended Newton–Raphson technique that is free of the above-described issues associated with the ordinary Newton–Raphson method. It is worthwhile to mention that the electrical analogy approach allows the setting of only one slack node for the entire MES model.

Finally, this paper is organized as follows: in Section 2, the developed electrical equivalent models of the most important natural gas network elements, as well as the electrical equivalent models of the linking facilities between the power and natural gas systems, are described. Two case studies are conducted in Section 3 to verify and demonstrate the applicability of the presented MES modeling approach, while conclusions are given in Section 4.

2. The Equivalent Electrical Analogy Models

In order to simulate the natural gas network, it is only important to know gas pressures, volumetric gas flows, and gas temperatures at the element suction (inlet) and discharge (outlet) side (if applicable). In other words, it is not relevant what happens inside the network element in detail, but what goes in/out of it. Hence, in this paper, the so-called network port theory and the load flow method formulation, known in the power systems analysis, were applied to define the electrical equivalent models of the gas network elements.

In the network port theory, the most used were the one- and two-port network models. The two-port network model, shown in Figure 1, has two-terminal pairs: input (port 1) and output (port 2), and it is represented by four external variables: input and output phase voltages (V_el1f and V_el2f, respectively), as well as input and output currents (I_el1 and I_el2, respectively) [34]. It should be noted that subscript _el in this paper was used to distinguish electrical and natural gas variables. Furthermore, the two-port network can be treated as a black box modeled by the relationships between the four external variables. Similarly, the one-port network has only an input port (one-terminal pair) and two external variables: V_el1f and I_el1 [35]. The currents flowing into the model were assumed to be positive, and the voltages had positive reference polarities [36].

In this paper, we differentiated two types of variables: internal and external. The internal variables referred to the element ports themselves, and the mandatory ones were, according to the load flow method [37], voltage magnitude V_el and angle A_el (in this case, the line voltage took place, so the index _f was omitted), as well as active P_el and reactive Q_el power. Additional variables may also be specified, such as the number of transformer taps, while some of them, such as the current flow, were derived. The external variables can refer to any other element in the model. Only the internal variables were mandatory. It should be noted that, as it has a strong influence on the form of equations, the model’s load flow equations together with their first partial derivatives by internal and external variables (elements of the Jacobian matrix) had to be specified in the load flow calculation process. Furthermore, it is common in the load flow formulation that the parameter values used in the load flow equations are expressed in the p.u. system (marked in bold) with the base power of 100 MVA. Moreover, these equations were expressed in an implicit form where the right part (=0) was omitted.

Since the reactive power does not exist in natural gas networks, the active power was equal to the apparent power, i.e., the power factor was equal to 1. Consequently, the three-phase active power can be calculated as:

$$P_{el}=\sqrt{3} · V_{el} · I_{el}$$

(1)

where V_el represents the line voltage.

It should be noted that in the p.u. system of the balanced power network, the phase and line voltages were equal [38,39].

The volumetric gas flow Q can be calculated as [32]:

$$Q=\frac{z·R_u·T}{M_w·p}·\dot{m}$$

(2)

where z, R_u, T, M_w, p, and $\dot{m}$ are the compressibility factor, the universal gas constant, the absolute gas temperature, the molecular weight of the natural gas composition, the gas pressure, and the gas mass flow, respectively.

The compressibility factor z can be calculated using several well-known equations. In this paper, Papay’s equation [40] was used:

$$z=1-3.52·p_r·\exp \left(-2.260·T_r\right)+0.274·p_r^2·\exp \left(-1.878·T_r\right)$$

(3)

where p_r and T_r are the reduced gas pressure and temperature, respectively.

Since the reduced gas pressure p_r was determined as a relation of the actual pressure to the critical value, to make the power flow equations shorter and easier to read, (3) was segmented into two constant parts. These parts were further, in the load flow equations, multiplied by the gas pressures (port voltages). Thus, for the two-port network element, the constant parts were:

$$k_{11}=3.52·k_{pr}·\exp \left(-2.26·T_1·k_{Tr}\right)$$

(4)

$$k_{21}=0.274·k_{pr}^2·\exp \left(-1.878·T_1·k_{Tr}\right)$$

(5)

$$k_{12}=3.52·k_{pr}·\exp \left(-2.26·T_2·k_{Tr}\right)$$

(6)

$$k_{22}=0.274·k_{pr}^2·\exp \left(-1.878·T_2·k_{Tr}\right)$$

(7)

where k₁₁ and k₂₁, as well as k₁₂ and k₂₂, are the first and the second constant parts of the first and second port compressibility factor equation, respectively. k_pr and k_Tr are constants calculated as the sum of all relations of gas portions in the natural gas mixture to its critical values (gas pressures and temperatures, respectively). T₁ and T₂ represent gas temperatures at both element sides.

In the following subsections, the load flow equations of the most important natural gas network elements, as well as P2G, and GFPP, are described.

2.1. Gas Compressor Station

The gas compressor station in this paper was modeled as the two-port network element. A detailed electrical analogy model of the gas compressor station can be found in [32].

In the load flow formulation, the gas compressor station was to be considered as the PV type, i.e., the active power and the voltage had to be specified/known.

Moreover, as the reactive power does not exist in the natural gas networks, the reactive power had to be set to zero.

The next set of equations depended on the compressor’s operation state, and they were based on the gas mass flow conservation through the compressor as well as on the gas pressures relation according to the compression ratio [32].

If the gas compressor was in operation, then the inlet gas mass flow minus the gas driver consumption had to be equal to the outlet gas mass flow, i.e.,:

$$\frac{p_2·Q_2}{z_2·T_2}=\frac{p_1·Q_1}{z_1·T_1}·\left(1-k_0/100\right)$$

(8)

where k₀ represents the gas mass flow consumption by the driver, expressed in % of the inlet gas mass flow.

By applying the electrical analogy approach to (8) and the described segmentation techniques to (3), and expressing the volumetric gas flow (i.e., electric current) using (1), the gas mass flow through the compressor station in p.u., i.e., the power equation, can be expressed as follows:

$$\begin{array}{c}{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{2}}·\left(1-k_{11}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kV{b}_1+k_{21}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}^{\mathbf{2}}·kV{b}_1^2\right)·T_1+{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}·(1-k_{12}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kV{b}_2+k_{22}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}^2·kV{b}_2^2)·T_2·\left(1-k_0/100\right)=0\hfill \end{array}$$

(9)

where kVb₁ and kVb₂ are inlet and outlet base pressures (voltages) in bar that should be applied to convert p.u. voltages into natural units as demanded by (3).

Moreover, the gas mass flow at the discharge side was flowing out of the two-port network model, so the positive sign came in front of the P_el1 part in (9).

If the compressor driver was electric-powered, then k₀ was equal to zero.

The discharge gas pressure, i.e., the voltage equation, depended on the compression ratio applied. In this paper, three types of compression ratios could be set. The first one was the user-specified fix compression ratio r_c, so the compressor station inlet and outlet gas pressures (voltages) related as:

$${\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kV{b}_2-r_c·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kV{b}_1=0$$

(10)

The second compression ratio option was to set the outlet pressure to the nominal value:

$${\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}-1.0=0$$

(11)

The third compression ratio option was to copy the p.u. pressure from the inlet to the outlet side:

$${\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}-{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}=0$$

(12)

If the gas compressor was bypassed, then the inlet and outlet sides had the same gas pressures, temperatures, and compressibility factors, while the driver consumption k₀ was equal to zero. The application of the former to (9) resulted in the inlet and outlet active powers equality (with opposite signs since discharge flow was out of the two-port network model):

$${\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{2}}+{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}=0$$

(13)

and in gas pressure equality (p.u. pressures had to be converted to the natural units by multiplying with their base values):

$${\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kV{b}_2-{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kV{b}_1=0$$

(14)

A gas compressor was also bypassed if it was electric-powered and the electric voltage of the bus from which the gas compressor driver was supplied (external variable) became an outage.

2.2. Compressor’s Electric Driver Consumption

If the gas compressor was electric-powered, then it was necessary to connect an appropriate load to the supplying node in the power network. In this paper, the electric load of the compressor driver was modeled as the one-port network element and as the PQ type in the load flow formulation (active and reactive powers had to be specified/known).

The electric load of the compressor driver needed to track the compressor operation state. Hence, the external variables needed to be defined as follows:

• V_el-ext1—compressor’s voltage at the suction side;
• P_el-ext1—compressor’s active power at the suction side;
• P_el-ext2—compressor’s active power at the discharge side.

In this paper, if the gas compressor was in operation, the reactive power was set to the fixed value:

$${\mathit{Q}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVAb-kq=0$$

(15)

where kq and kVAb are the reactive power constant settings by the user in kvar and base power of 10⁵ kVA, respectively.

Moreover, the driver’s active power consumption could be determined in two ways. The first one was to calculate the active power consumption as a percentage of the gas mass flow through the compressor. Accordingly, to determine the power consumption, gas energy had to be converted into kWh using the gas heating value. In some countries, a low heating value (LHV) is used, while in others, a high heating value (HHV) takes place. So, in this paper, the heating value was defined as a constant that could be set by the user. The volumetric gas flow at the compressor’s suction side could be determined using the defined external variables V_el-ext1, P_el-ext1, and (1). Since the gas heating value was given per standard gas volume, the calculated actual volumetric gas flow in A, i.e., in m³/s, needed to be converted to the standard gas flow (gas flow at standard conditions: 15 °C and 101.325 kPa, in Nm³/s) as follows:

$$Q_{st}=Q_1·\frac{p_1}{p_{st}}·\frac{z_{st}}{z_1}·\frac{T_{st}}{T_1}$$

(16)

where subscripts _st and ₁ denote standard conditions and the compressor’s station suction side, respectively.

The suction gas pressure p₁ was represented with V_el-ext1, and the compressibility factor z₁ was calculated using the segmentation of (3), i.e., multiplying (4) and (5) with V_el-ext1. Finally, the compressor driver power consumption P_el1 was defined as:

$$\begin{array}{c}{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}-\frac{{\mathit{P}}_{\mathit{e}\mathit{l}\text{-}\mathit{e}\mathit{x}\mathit{t}\mathbf{1}}}{\sqrt{3}}·\frac{{10}^5}{p_{st}}·\frac{z_{st}}{\left(1-k_{11}·{\mathit{V}}_{\mathit{e}\mathit{l}\text{-}\mathit{e}\mathit{x}\mathit{t}\mathbf{1}}·kVb+k_{21}·{\mathit{V}}_{\mathit{e}\mathit{l}\text{-}\mathit{e}\mathit{x}\mathit{t}\mathbf{1}}^{\mathbf{2}}·kV{b}^2\right)}·\frac{T_{st}}{T_1}·HV·3600\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}·k_0/100=0\hfill \end{array}$$

(17)

It should be noted that the gas pressure p₁ in (16) was in Pa, while in the load flow equation the p.u. voltage multiplied by the base value gave pressure in bar. Because of that, factor 10⁵ exists in (17). Moreover, factor 3600 in (17) was a result of the conversion of kWh/h into kWh/s.

The second possibility was to set the active power manually to the fixed value:

$${\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVAb-kp=0$$

(18)

where kp is the active power constant setting by the user in kW.

If the gas compressor station inlet and outlet active powers in the unique electrical MES model (P_el-ext1 and P_el-ext2) were equal, then the compressor was bypassed, and, consequently, the driver’s electric power consumption (P_el1 and Q_el1) equaled zero.

2.3. Gas Pressure Reducing and Metering Station

The gas pressure reducing and metering station (GPRMS) in this paper was also modeled (such as the gas compressor station) as the two-port network element and the PV type in the load flow formulation. A detailed GPRMS electrical analogy model can be found in [41].

As for the gas compressor, the GPRMS reactive power also had to be set to zero.

Regarding the gas mass flow balance through the GPRMS, two possibilities exist: GPRMS with or without gas preheating prior to entering the control valve [41]. In the first case, when the preheating process took place, the gas mass flow in p.u., i.e., the power equation, through GPRMS was similar as for the gas compressor, but without the gas compressor driver’s consumption, i.e.,:

$$\begin{array}{c}{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{2}}·\left(1-k_{11}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kV{b}_1+k_{21}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}^{\mathbf{2}}·kV{b}_1^2\right)·T_1+{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}·(1-k_{12}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kV{b}_2+k_{22}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}^{\mathbf{2}}·kV{b}_2^2)·T_2=0\hfill \end{array}$$

(19)

In the second case, when the gas was not preheated prior to entering the control valve, the Joule–Thomson effect had to be taken into account. Accordingly, the discharge gas temperature can be calculated as follows [41]:

$$T_2=T_1 - \left({\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kV{b}_1-{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}} · kV{b}_2\right)·0.08058·JT$$

(20)

where JT is the Joule–Thomson coefficient in °F/100 psi.

In this paper, JT of 7 °F/100 psi, as suggested by the American Gas Association (AGA), was used, but the user could set any value to the JT variable in the model.

Since the discharge gas temperature Equation (20) contains suction and discharge gas pressures, the discharge compressibility factor equation had to be segmented into six constant terms k_1v2–k_6v2:

$$k_{1v2}=3.52·k_{pr}·\exp \left(-2.26·k_{Tr}·T_1\right)$$

(21)

$$k_{2v2}=0.182103·k_{Tr}·JT$$

(22)

$$k_{3v2}=0.274·k_{pr}^2·\exp \left(-1.878·k_{Tr}·T_1\right)$$

(23)

$$k_{4v2}=0.151323·k_{Tr}·JT$$

(24)

$$k_{5v2}=T_1$$

(25)

$$k_{6v2}=0.08058·JT$$

(26)

Applying (20)–(26) into (19), the gas mass balance through GPRMS without a preheating process resulted in:

$$\begin{array}{c}{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{2}}·\left(1-k_{11}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kV{b}_1+k_{21}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}^{\mathbf{2}}·kV{b}_1^2\right)·T_1+{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}·(1-k_{1v2}·\exp (k_{2v2}·({\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}· kV{b}_1-{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kV{b}_2\left)\right)·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kV{b}_2+k_{3v2}·\exp (k_{4v2}·({\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kV{b}_1\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}-{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kV{b}_2\left)\right)·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}^{\mathbf{2}}·kV{b}_2^2)·(k_{5v2}-k_{6v2}·({\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kV{b}_1-{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}· kV{b}_2\left)\right)=0\hfill \end{array}$$

(27)

In both cases, the discharge pressure was set to the nominal value, so (11) could be applied as the voltage equation.

If the suction gas pressure was lower than the nominal discharge value, then the GPRMS was bypassed, so (13) and (14) could be applied.

2.4. Gas Pipeline

The electrical analogy model of the gas pipeline in a steady state is given in [33], and it is represented by a serial connection of electrical resistance and the constant electric potential source. The former represents the resistance to the gas flow through the pipeline, while the latter imports the effect of the pipeline inclination, i.e., potential energy into the electrical analogy. If the pipeline is horizontal, the latter part of the model does not exist.

However, in [33], the standard volumetric gas flow is analog to the electric current, while in this paper, the actual volumetric gas flow took place since it was directly proportional to the gas velocity upon which the gas flow was constrained through the pipeline. Thus, it was necessary to derive the resistance equation in terms of the actual gas flow and p.u. values.

The gas pressure drop along the pipeline can be calculated as [33]:

$$p_1^2-p_2^2-E_p=R_p·Q_{st}$$

(28)

where E_p and R_p are the constant electrical potential source, i.e., potential energy due to the pipeline inclination, and pipeline resistance, respectively.

The potential energy E_p can be determined as follows [33]:

$$E_p=0.06843·G · \left(H_2-H_1\right) · \frac{p_{avg}^2}{T_{avg} · z_{avg}}$$

(29)

where H₁ and H₂ are heights of the pipe upstream and downstream side, respectively. Subscript _avg represents the average value between the upstream and downstream variables.

Moreover, in [33], several equations for the gas pipeline resistance in a fully turbulent flow are given, according to the Darcy friction factor expression applied. In this paper, the AGA fully turbulent gas pipeline resistance equation was used, as a tradeoff between the accuracy of the calculation results and the simplicity of expression:

$$R_p=\frac{p_{st}^2·L·G·T_{avg}·z_{avg}·Q_{st}}{176.85·T_{st}^2·{\eta }_p^2·\left(p_1+p_2\right)·D^5}·{\left[-2·{\log }_{10}\left(\frac{\epsilon }{3.7·D}\right)\right]}^{-2}$$

(30)

where L, G, η_p, D, and ε are the pipeline length, gas specific gravity, pipeline efficiency factor, pipeline diameter, and surface roughness, respectively.

The gas specific gravity G can be substituted with molecular weight M_w as follows [33]:

$$G=\frac{M_w}{M_{w air}}\approx \frac{M_w}{29}$$

(31)

The pipeline efficiency factor η_p represents additional friction imposed by fittings (e.g., valves, bends, tees, etc.), corrosion, dust deposition, etc. [33].

To convert the standard to actual volumetric gas flow (28), (30) and (16) should be applied.

The gas pipeline in this paper was, such as the compressor station and GPRMS, modeled as the two-port network element and the PV-type in the load flow formulation. The reactive power of both pipeline sides also has to be set to zero.

After applying the segmentation to the compressibility factors using (4)–(7) and the substitution of the current with active power and voltage using (1) in (30), and applying (30) to (28), the pressure equation is as follows:

$$\begin{array}{c}{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}^{\mathbf{2}}-{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}^{\mathbf{2}}-E_{p const}·{\left({\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}+{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}\right)}^2/(2-k_{11}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVb+k_{21}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}^{\mathbf{2}}·kV{b}^2-k_{12}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kVb+k_{22}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}^{\mathbf{2}}·kV{b}^2)-{\mathit{R}}_{\mathit{p} \mathit{c}\mathit{o}\mathit{n}\mathit{s}\mathit{t}}/6·(2-k_{11}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVb\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}+k_{21}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}^{\mathbf{2}}·kV{b}^2-k_{12}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kVb+k_{22}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}^{\mathbf{2}}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}·kV{b}^2)/{\left(1-k_{11}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVb+k_{21}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}^{\mathbf{2}}·kV{b}^2\right)}^2·\left|{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}\right|·{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVAb\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}·{10}^5=0\hfill \end{array}$$

(32)

where R_{p const} and E_{p const} represent constants calculated as follows:

$${\mathit{R}}_{\mathit{p} \mathit{c}\mathit{o}\mathit{n}\mathit{s}\mathit{t}}=\frac{L·M_w·\left(T_1+T_2\right)·z_{st}^2}{10257.46·{\mathsf{\eta }}_p^2·D^5·T_1^2}·\frac{1}{{\left(-2·{\log }_{10}\left(\frac{\mathsf{\epsilon }}{3.7·D}\right)\right)}^2}·\frac{kVAb}{kV{b}^2}·{10}^5$$

(33)

$$E_{p const}=\frac{M_w · \left(H_2-H_1\right)}{423.79 · \left(T_1+T_2\right)}$$

(34)

Moreover, the absolute value of P_el1 in (32) was needed as the gas pipeline model had to be reciprocal, i.e., upstream and downstream sides could be freely switched in the network model. The same situation could also be observed when the gas flow changed its direction through the pipeline.

The gas mass flow had to be equal at both pipeline sides, i.e., the mass flow equality must be satisfied. Hence, applying the segmentation of z and (1) into (2), the power equation is as follows:

$$\begin{array}{c}{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}·\left(1-k_{12}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kVb+k_{22}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}^{\mathbf{2}}·kV{b}^2\right)·T_2+{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{2}}·(1-k_{11}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVb+k_{21}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}^{\mathbf{2}}·kV{b}^2)·T_1=0\hfill \end{array}$$

(35)

2.5. Gas Load, Gas Input, Underground Gas Storage

In this paper, gas load, gas input, and underground gas storage were modeled as the one-port network elements and the PQ type in the load flow formulation.

The below-derived load flow models were similar for all three gas network elements. The difference was only in the sign of the desired volumetric gas flow that is usually defined at standard conditions Q_st and given in 1000∙Nm³/h. For the gas load, a positive value was used, for the gas input a negative value was used, while for the underground gas storage both signs were possible, depending on the operation state (positive for the gas injection process and negative for the gas withdrawal process).

Having defined Q_st, the gas mass flow $\dot{m}$ in kg/s can be calculated as follows [32]:

$$\dot{m}=\frac{Q_{st}·p_{st}·Mw}{3600·z_{st}·R_u·T_{st}}$$

(36)

As in all the natural gas network elements described above, the reactive power had to be set to zero.

The power equation, derived combining (1) and (2), is as follows:

$${\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}-\sqrt{3}·\left(1-k_{11}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVb+k_{21}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}^{\mathbf{2}}·kV{b}^2\right)·\frac{\dot{m}·R_u·T_1}{M_w}·\frac{1}{100·kVAb}=0$$

(37)

2.6. Power-to-Gas

The power-to-gas facility in this paper was also modeled as the two-port network element and the PQ type in the load flow formulation. The first port was dedicated to the power network side, while the second port was connected to the natural gas network side. The input values were power consumption kp in MW, the overall process efficiency η_P2G in %, as well as the temperature of the produced gas T₂ in K.

Since the P2G facility consumed only active power and bearing in mind that reactive power does not exist in natural gas networks, the reactive power of both ports had to be set to zero.

In the power equation of the first port, the active power was set to the defined fixed value:

$${\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVAb-kp·1000=0$$

(38)

In the power equation of the second port, the active power was calculated using (1), (3), (6), (7) and (16) as follows:

$${\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{2}}+\sqrt{3}·\frac{{\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}·p_{st}·T_2}{3600·HV·z_{st}·T_{st}}·\left(1-k_{12}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}·kVb+k_{22}·{\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{2}}^{\mathbf{2}}·kV{b}^2\right)·\frac{{\eta }_{P2G}}{{10}^7}=0$$

(39)

2.7. Gas-Fired Power Plant

As power-to-gas, the gas-fired power plant (GFPP) was also modeled as the two-port network element where the first node was connected to the power network, while the second node was connected to the natural gas network.

The gas temperature T₂ in K, the overall power plant efficiency η_GFPP in %, and the power generation kp in MW represented the input values.

Furthermore, there were two possibilities in the GFPP operation regime. In the first one, the voltage of the power bus at which the GFPP was connected could be regulated, while in the second one, a fixed value of the reactive power generated kq in Mvar could be set. In the load flow formulation, the former was represented as the PV type, while the latter was represented as the PQ type.

Moreover, the GFPP model defined the minimum allowed gas pressure of the gas node at which GFPP was connected. If the threshold setting was violated, the GFPP was disconnected, so the active and reactive powers of both GFPP’s sides had to be set to zero.

If the GFPP regulated the voltage of the connecting bus to the desired value kv expressed in % of the nominal value, the voltage equation was applied:

$${\mathit{V}}_{\mathit{e}\mathit{l}\mathbf{1}}-\frac{kv}{100}=0$$

(40)

Moreover, if the fix reactive power generation kq in Mvar was set, the reactive power equation was executed:

$${\mathit{Q}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVAb+kq·1000=0$$

(41)

Since the reactive power does not exist in the natural gas network, it had to be set to zero at the second port.

The active power equation of the first port was similar as for power-to-gas, but with the opposite sign:

$${\mathit{P}}_{\mathit{e}\mathit{l}\mathbf{1}}·kVAb+kp·1000=0$$

(42)

Finally, the active power equation of the second port was the same as (39), but the variable η_GFPP instead of η_P2G should be used.

3. Case Studies

In this section, two case studies were conducted to verify the presented equivalent electrical analogy models and applied modeling approach described in the previous section.

The presented equivalent electrical analogy models were written in the C/C++ scripts and compiled to the Dynamic Link Library (DLL) files, which could be further assigned in the NEPLAN (a well-known software package for power network simulations) graphical editor as the User-Defined Models (UDM) [42]. The C/C++ user’s code must contain several routines and follow certain instructions that are described in the NEPLAN manual [42].

In the first case study, an extensive gas network was modeled in the NEPLAN electric module, using the presented equivalent electrical models of the gas network elements, as well as in the well-known commercial software for the natural gas network analyses—SIMONE. The verification of the results obtained by the NEPLAN was conducted by a comparison with the results obtained by SIMONE.

After the verification of the accuracy of the applied methodology in the first case study, in the second case study, an electrical power network was added to the same model together with the main coupling facilities between the power and gas networks: the GFPP, the P2G, and the electric-driven natural gas compressor station. The goal of the second case study was to demonstrate the applicability of the presented modeling approach as well as to emphasize the importance of having the unique MES simulation model.

The same gas mixture from

Table 1. Composition of the processed natural gas.

Gas	Formula	Mole Fraction	Molecular Weight (g/mol)	Molecular Weight of Gas Component (g/mol)
Methane	CH4	95.2%	16.043	15.34
Ethane	C2H6	2.29%	30.07	0.69
Propane	C3H8	0.63%	44.096	0.28
i-butane	C4H10	0.09%	58.124	0.05
Nitrogen	N2	1.06%	28.014	0.30
Carbon-dioxide	CO2	0.31%	44.01	0.14
Molecular weight of the gas mixture				16.79

was used in all simulations within both case studies. Additionally, in all calculations, the gas heating value of 9.260625 kWh/Nm³ was used. Reference conditions were 15 °C and 101.325 kPa.

3.1. Case Study 1

The gas network used to verify the applied approach represented an imaginary network consisting of 35 nodes, 27 pipelines, 4 gas compressor stations (CS), 6 GPRMSs, 7 gas exits/loads, 2 gas inputs, and 1 underground gas storage (UGS). It was designed to cover a wide range of input data: gas inputs and loads, pipeline diameters, roughness, lengths, and elevations. The model of the gas network, constructed in SIMONE, is shown in Figure 2. Two valves VA1 and VA2 served to switch off the pipelines N2–N7 and N7–N8 out of operation in situations when the gas flow velocities through them fell below the critical value. In those situations, all downstream natural gas flows through the pipeline N2–N8.

Within case study one, six simulations were performed, as indicated in

Table 2. Scenario definition in case study 1.

UGS	CS1	CS2	Scenario
Injection	In operation	In operation	S1
	Bypass	In operation	S2
	In operation	Bypass	S3
Withdrawal	In operation	In operation	S4
	Bypass	In operation	S5
	In operation	Bypass	S6

. The same gas pipeline parameters, presented in

Table 3. Gas pipeline parameters.

Node 1	Node 2	Diameter (mm)	Roughness (mm)	Length (km)	Node 1 Height (m)	Node 2 Height (m)
N18	EXIT4	200	0.018	10	100	120
INPUT1	N1	600	0.012	10	0	10
N1	N2	600	0.012	2	10	10
N17	N18	200	0.018	5	4	100
N2	N8	600	0.012	70	10	10
N13	N14	150	0.012	30	80	80
N18	EXIT5	200	0.018	5	100	80
N14	EXIT3	150	0.012	5	80	150
N20	N19	200	0.03	50	25	50
N8	EXIT2	400	0.012	50	10	15
INPUT2	N20	200	0.012	4	0	25
N21	N22	400	0.012	20	70	100
EXIT2	N9	400	0.012	50	15	0
N9	N12	150	0.012	100	0	80
N23	EXIT7	350	0.012	10	100	100
N10	N11	400	0.012	10	0	0
N2	N7	350	0.015	15	10	−10
N7	N8	350	0.015	20	−10	10
N15	N19	400	0.018	4	30	50
N2	N3	500	0.012	80	10	20
N23	N24	250	0.03	20	100	65
N15	N16	200	0.012	4.5	30	4
N8	N15	400	0.012	10	10	30
N5	N6	350	0.015	25	100	100
N24	N25	250	0.01	20	65	0
N4	N5	350	0.024	32	20	100
N19	N21	400	0.012	20	50	70

, were used in all defined scenarios. The pipeline efficiency was taken as 1, and the total length of the gas network was 681.5 km.

The simulation parameters, shown in

Table 4. Simulation parameters.

Element	Parameter	UGS Injection	UGS Withdrawal	Unit
CS1	Discharge pressure	75	75	bar
CS1	Discharge temperature	15	15	°C
CS2	Discharge pressure	75	75	bar
CS2	Discharge temperature	20	20	°C
CS3	Compression ratio	Off	1.5	-
CS3	Discharge temperature	10	10	°C
CS4	Discharge pressure	125	Off	bar
CS4	Discharge temperature	10	10	°C
EXIT1	Gas demand	200	200	1000 Nm3/h
EXIT2	Gas demand	35	35	1000 Nm3/h
EXIT3	Gas demand	3	3	1000 Nm3/h
EXIT4	Gas demand	20	20	1000 Nm3/h
EXIT5	Gas demand	15	15	1000 Nm3/h
EXIT6	Gas demand	77	77	1000 Nm3/h
EXIT7	Gas demand	155	155	1000 Nm3/h
GPRMS1	Discharge pressure	35	35	bar
GPRMS1	Discharge temperature	10	10	°C
GPRMS2	Discharge pressure	35	35	bar
GPRMS2	Discharge temperature	10	10	°C
GPRMS3	Discharge pressure	35	35	bar
GPRMS3	Discharge temperature	10	10	°C
GPRMS4	Discharge pressure	30	30	bar
GPRMS4	Discharge temperature	18	18	°C
GPRMS5	Discharge pressure	30	Off	bar
GPRMS5	Discharge temperature	10	10	°C
GPRMS6	Suction pressure	Off	125	bar
GPRMS6	Discharge temperature	10	10	°C
INPUT1	Pressure set	75	75	Bar
INPUT1	Gas temperature	10	10	°C
INPUT2	Gas supply	35	−15	1000 Nm3/h
INPUT2	Gas temperature	10	10	°C
UGS	Gas demand	80	−150	1000 Nm3/h

, differ for the UGS injection and the withdrawal process. When the UGS is in the gas injection phase, the gas reduction station GPRMS5 and the compressor station CS4 were in operation. Otherwise, when natural gas was withdrawn from the UGS, the GPRMS6 and the CS3 were in operation. During the withdrawing process, the second gas source (i.e., connection to another gas network) INPUT2 changed the role and became a gas consumer. For all compressors, the driver’s gas consumptions were neglected; therefore, in this case, they could be considered as electric-driven compressors.

The main goal was to determine node gas pressures and branch volumetric gas flows. The deviations between the obtained results were calculated as:

$$\mathrm{Deviation}=\frac{\mathrm{NEPLAN}-\mathrm{SIMONE}}{\mathrm{SIMONE}}·100\%$$

(43)

where NEPLAN and SIMONE denote results obtained by models developed in NEPLAN and SIMONE, respectively.

The results of the node gas pressures and branch volumetric gas flows are shown in

Table 5. Results of simulation in case study 1—gas pressures (bar).

Node	S1-N	S1-S	S2-N	S2-S	S3-N	S3-S	S4-N	S4-S	S5-N	S5-S	S6-N	S6-S
EXIT1	35.000	35.000	28.126	27.200	35.000	35.000	35.000	35.000	31.300	30.556	35.000	35.000
EXIT2	66.286	65.602	66.236	65.604	66.218	65.629	75.728	75.697	75.699	75.698	75.678	75.714
EXIT3	34.259	34.125	34.259	34.125	34.259	34.125	34.259	34.125	34.259	34.125	34.259	34.125
EXIT4	30.044	29.795	30.044	29.795	30.044	29.795	30.044	29.795	30.044	29.795	30.044	29.795
EXIT5	31.544	31.386	31.544	31.386	31.544	31.386	31.544	31.386	31.544	31.386	31.544	31.386
EXIT6	30.000	30.000	30.000	30.000	19.630	20.213	30.000	30.000	30.000	30.000	23.542	24.539
EXIT7	72.652	72.564	72.652	72.564	48.660	47.934	72.652	72.564	72.652	72.564	50.264	49.747
INPUT1	75.000	74.995	75.000	74.995	75.000	74.995	75.000	74.996	75.000	74.996	75.000	74.996
INPUT2	79.935	79.818	79.895	79.820	79.808	79.867	64.693	64.140	64.658	64.142	64.519	64.200
N1	73.173	73.095	73.137	73.097	73.157	73.101	74.145	74.098	74.120	74.099	74.134	74.102
N10	30.000	30.000	30.000	30.000	30.000	30.000	53.764	53.868	53.746	53.868	53.733	53.878
N11	29.163	29.053	29.163	29.053	29.163	29.053	55.271	55.445	55.253	55.446	55.241	55.455
N12	63.256	62.108	63.204	62.110	63.185	62.136	79.478	79.484	79.451	79.485	79.432	79.500
N13	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000
N14	34.524	34.410	34.524	34.410	34.524	34.410	34.524	34.410	34.524	34.410	34.524	34.410
N15	67.042	66.579	66.994	66.581	66.912	66.629	68.701	68.369	68.669	68.370	68.569	68.414
N16	66.118	65.589	66.069	65.591	65.985	65.641	67.811	67.418	67.778	67.419	67.676	67.464
N17	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000	35.000
N18	32.026	31.914	32.026	31.914	32.026	31.914	32.026	31.914	32.026	31.914	32.026	31.914
N19	66.041	65.538	65.991	65.540	65.885	65.598	67.220	66.841	67.187	66.842	67.055	66.897
N2	72.814	72.721	72.770	72.723	72.794	72.729	73.985	73.930	73.956	73.931	73.971	73.934
N20	78.992	78.833	78.951	78.834	78.864	78.882	64.740	64.209	64.706	64.211	64.568	64.268
N21	59.797	58.963	59.742	58.965	59.478	59.087	61.108	60.427	61.071	60.429	60.779	60.539
N22	52.696	51.417	52.633	51.420	52.163	51.625	54.193	53.114	54.151	53.115	53.652	53.299
N23	75.000	75.000	75.000	75.000	52.163	51.623	75.000	75.000	75.000	75.000	53.652	53.297
N24	67.012	66.875	67.012	66.875	38.193	38.234	67.012	66.875	67.012	66.875	40.247	40.549
N25	59.969	59.454	59.969	59.454	19.630	20.214	59.969	59.454	59.969	59.454	23.542	24.540
N3	67.669	67.259	67.317	67.274	67.648	67.267	68.937	68.578	68.606	68.590	68.922	68.583
N4	75.000	75.000	67.317	67.272	75.000	75.000	75.000	75.000	68.606	68.587	75.000	75.000
N5	58.981	58.715	48.612	48.521	58.981	58.715	58.981	58.715	50.423	50.376	58.981	58.715
N6	44.199	43.392	28.126	27.201	44.199	43.390	44.199	43.392	31.300	30.557	44.199	43.390
N7	71.692	71.528	71.647	71.530	71.653	71.543	0.000	0.000	0.000	0.000	0.000	0.000
N8	69.945	69.615	69.899	69.617	69.881	69.640	72.792	72.608	72.762	72.609	72.741	72.626
N9	64.485	63.540	64.434	63.543	64.414	63.568	80.645	80.801	80.618	80.802	80.599	80.817
UGS	125.000	125.000	125.000	125.000	125.000	125.000	125.000	125.000	125.000	125.000	125.000	125.000

N—NEPLAN; S—SIMONE.

and

Table 6. Results of simulation in case study 1—volumetric gas flows (m³/s).

Node	S1-N	S1-S	S2-N	S2-S	S3-N	S3-S	S4-N	S4-S	S5-N	S5-S	S6-N	S6-S
EXIT1	1.448	1.448	1.832	1.898	1.448	1.448	1.448	1.448	1.633	1.658	1.448	1.448
EXIT2	0.296	0.306	0.297	0.306	0.297	0.305	0.451	0.449	0.452	0.449	0.452	0.449
EXIT3	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022
EXIT4	0.171	0.172	0.171	0.172	0.171	0.172	0.171	0.172	0.171	0.172	0.171	0.172
EXIT5	0.121	0.122	0.121	0.122	0.121	0.122	0.121	0.122	0.121	0.122	0.121	0.122
EXIT6	0.682	0.682	0.682	0.682	1.067	1.029	0.682	0.682	0.682	0.682	0.881	0.844
EXIT7	0.499	0.500	0.499	0.500	0.782	0.795	0.499	0.500	0.499	0.500	0.754	0.763
INPUT1	1.707	1.716	1.725	1.715	1.715	1.713	1.148	1.154	1.166	1.153	1.157	1.151
INPUT2	0.101	0.101	0.101	0.101	0.101	0.101	0.055	0.056	0.055	0.056	0.055	0.056
N1	1.756	1.767	1.775	1.766	1.764	1.763	1.164	1.170	1.182	1.169	1.172	1.167
N10	0.684	0.697	0.684	0.697	0.684	0.697	0.677	0.676	0.678	0.676	0.678	0.676
N11	0.705	0.722	0.705	0.722	0.705	0.722	0.657	0.655	0.657	0.655	0.657	0.654
N12	0.011	0.012	0.011	0.012	0.011	0.012	0.009	0.009	0.009	0.009	0.009	0.009
N13	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022
N14	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022	0.022
N15 (N16)	0.123	0.124	0.123	0.124	0.124	0.124	0.120	0.121	0.120	0.121	0.120	0.121
N15 (N19)	0.694	0.703	0.695	0.703	0.705	0.699	0.847	0.855	0.847	0.855	0.858	0.851
N15 (N8)	0.818	0.828	0.818	0.828	0.829	0.823	0.967	0.975	0.967	0.975	0.978	0.971
N16	0.125	0.126	0.125	0.126	0.126	0.126	0.122	0.123	0.122	0.123	0.122	0.122
N17	0.253	0.253	0.253	0.253	0.253	0.253	0.253	0.253	0.253	0.253	0.253	0.253
N18 (EXIT4)	0.159	0.160	0.159	0.160	0.159	0.160	0.159	0.160	0.159	0.160	0.159	0.160
N18 (EXIT5)	0.120	0.120	0.120	0.120	0.120	0.120	0.120	0.120	0.120	0.120	0.120	0.120
N18 (N17)	0.279	0.280	0.279	0.280	0.279	0.280	0.279	0.280	0.279	0.280	0.279	0.280
N19 (N15)	0.706	0.716	0.707	0.716	0.717	0.712	0.868	0.877	0.869	0.877	0.880	0.873
N19 (N20)	0.125	0.127	0.126	0.127	0.126	0.126	0.053	0.053	0.053	0.053	0.053	0.053
N19 (N21)	0.832	0.842	0.832	0.842	0.843	0.838	0.815	0.824	0.816	0.824	0.827	0.820
N2 (N1)	1.766	1.777	1.785	1.776	1.774	1.774	1.166	1.173	1.185	1.172	1.175	1.170
N2 (N3)	0.642	0.644	0.661	0.643	0.642	0.644	0.630	0.632	0.649	0.631	0.631	0.632
N2 (N7)	0.290	0.292	0.291	0.292	0.293	0.291	0.000	0.000	0.000	0.000	0.000	0.000
N2 (N8)	0.833	0.842	0.834	0.842	0.839	0.839	0.536	0.541	0.536	0.541	0.544	0.538
N20	0.102	0.103	0.102	0.103	0.103	0.103	0.055	0.056	0.055	0.056	0.055	0.055
N21	0.930	0.949	0.931	0.949	0.946	0.943	0.908	0.923	0.908	0.923	0.924	0.917
N22	1.071	1.106	1.073	1.106	1.096	1.096	1.039	1.066	1.039	1.066	1.062	1.058
N23 (EXIT7)	0.509	0.509	0.509	0.509	0.724	0.758	0.509	0.492	0.509	0.492	0.702	0.707
N23 (N24)	0.253	0.253	0.253	0.253	0.372	0.377	0.253	0.253	0.253	0.253	0.360	0.364
N24	0.272	0.272	0.272	0.272	0.524	0.506	0.272	0.272	0.272	0.272	0.495	0.475
N25	0.308	0.311	0.308	0.311	1.067	1.001	0.308	0.311	0.308	0.311	0.881	0.816
N3	0.698	0.703	0.721	0.702	0.698	0.703	0.683	0.688	0.706	0.687	0.683	0.688
N4	0.639	0.639	0.721	0.715	0.639	0.639	0.639	0.639	0.706	0.699	0.639	0.639
N5	0.814	0.818	1.010	1.012	0.814	0.819	0.814	0.818	0.970	0.954	0.814	0.819
N6	1.122	1.145	1.832	1.898	1.122	1.145	1.122	1.145	1.633	1.626	1.122	1.145
N7	0.296	0.297	0.296	0.297	0.298	0.296	0.000	0.000	0.000	0.000	0.000	0.000
N8 (EXIT2)	0.397	0.404	0.397	0.404	0.397	0.404	0.360	0.359	0.360	0.359	0.360	0.358
N8 (N15)	0.779	0.787	0.780	0.787	0.789	0.783	0.906	0.911	0.906	0.911	0.915	0.908
N8 (N2)	0.872	0.884	0.873	0.884	0.879	0.882	0.546	0.552	0.546	0.552	0.555	0.549
N8 (N7)	0.304	0.307	0.304	0.307	0.307	0.306	0.000	0.000	0.000	0.000	0.000	0.000
N9 (EXIT2)	0.306	0.317	0.306	0.317	0.306	0.317	0.420	0.417	0.420	0.417	0.420	0.417
N9 (N12)	0.011	0.011	0.011	0.011	0.011	0.011	0.009	0.009	0.009	0.009	0.009	0.009
UGS	0.139	0.139	0.139	0.139	0.139	0.139	0.261	0.261	0.261	0.261	0.261	0.261

N—NEPLAN; S—SIMONE.

, respectively. It should be noted that for each node, Kirchhoff’s current law had to be satisfied. This meant that the sum of all volumetric gas flows going into a node had to be equal to the sum of all volumetric gas flows going out of the same node. Since the gas temperature was taken as a constant along the pipeline and the pipe cross-section remained unchanged, according to (2) the volumetric gas flow at the pipe exit was greater than at the pipe inlet side, as a result of the gas pressure drop. For this reason, both the pipe inlet and the outlet volumetric gas flows were compared to the SIMONE results. If some node only had two connected branches, then the volumetric gas flows of both branches at the sides connected to that node were the same. In that case, only one value was presented in Table 6. However, if some node had three or more connected branches, then Table 6 contains the volumetric gas flows of all connected branches. In those situations, the name of the second node of the corresponding branch was specified in parentheses in order to provide information about the volumetric gas flow branch affiliation.

Due to different equations used for the calculation of the gas pressure drop along the gas pipeline in SIMONE and the proposed equivalent electrical model in NEPLAN, it was evident that deviations between the obtained results would exist. The node gas pressure deviations, as well as the branch volumetric gas flow deviations, were calculated and presented in Figure 3. In normal situations (scenarios S1 and S4), the absolute node gas pressure deviations were up to 2.49% (standard deviation 0.56%), while the absolute branch volumetric gas flow deviations were up to 3.50% (standard deviation 0.82%). Larger deviations were recorded in the N-1 situations when CS1 or CS2 were out of operation, i.e., bypassed. However, the standard deviations were slightly higher and were still very acceptable (0.81% and 1.18% for the node gas pressure and branch volumetric gas flow deviations, respectively). Hence, it can be concluded that the developed equivalent electrical models could efficiently simulate natural gas network steady-state load flows.

3.2. Case Study 2

In this case study, the electrical network, together with the coupling facilities, was added, so a unique electrical model of the MES, as shown in Figure 4, was constructed. The power network consisted of 14 nodes, 13 loads, 12 lines, 2 GFPPs, 1 WPP, and 2 transformers that couples 2 voltage levels of 110 kV and 20 kV. The coupling facilities between the gas and the power network were two GFPPs, the P2G and the electric driven gas compressor CS2, marked with bold lines in Figure 4. It should be noted that the GFPP1el and the GFPP2el were numerically equal to the GFPP1 and the GFPP2, respectively, and should have been used only when any of the networks were going to be analyzed separately, such as in case study one. In those situations, GFPP1 and GFPP2 should be disconnected. The frequency of the power system was taken as 50 Hz.

The parameters of electric lines and loads were given in

Table 7. Electric line parameters.

Node 1	Node 2	Resistance (Ω)	Reactance (Ω)	Susceptance (μS)	Conductivity (μS)	Rated Current (A)	Length (km)
B1	B2	0.164	0.338	3.456	0	450	50
B1	B3	0.164	0.338	3.456	0	450	50
B3	B6	0.164	0.338	3.456	0	450	50
B4	B5	0.164	0.338	3.456	0	450	50
B5	B3	0.164	0.338	3.456	0	450	50
B8	B7	0.3	0.41	2.011	0	290	5
B8	B9	0.3	0.41	2.011	0	290	5
B9	B10	0.3	0.41	2.011	0	290	5
B9	B11	0.3	0.41	2.011	0	290	5
B12	B11	0.3	0.41	2.011	0	290	5
B12	B14	0.3	0.41	2.011	0	290	5
B13	B12	0.3	0.41	2.011	0	290	5

and

Table 8. Electric loads.

Name	Active Power (MW)	Reactive Power (Mvar)
L2	20	2
L3	10	2
L4	16	3
L5	15	2
L6	60	5
L7	1.5	0.2
L8	3	1
L9	6	1
L10	12	2.5
L11	3	1
L12	6	2
L13	15	3
L14	3.5	1

, respectively. The total length of the power network was 285 km, and all lines were overhead types.

Each of the two identical 110/20 kV transformers, with parameters presented in

Table 9. Transformer parameters.

Name	Vector Group	Rated Power (MW)	Rated Primary Voltage (kV)	Rated Secondary Voltage (kV)	Short-Circuit Voltage (%)	Copper Losses (%)	Open-Circuit Current (%)	Iron Core Losses (kW)
TRF1	YNd5	40	110	20	17.59	0.25	0.51	20.2
TRF2	YNd5	40	110	20	17.59	0.25	0.51	20.2
TR-slack	YNd5	80	110	75	17.59	0	0	0

, on the primary side had an on-load tap changer with a total of 21 steps of 1.5% of nominal voltage. The secondary voltage was regulated to 103% of the nominal voltage. Moreover, one of the advantages of a unique electrical MES model was the opportunity to have only one slack node for both networks. For this reason, the third (imaginary) transformer, called TR-slack, was added into the unique model to couple busbars INPUT1 and B1 having the network feeders (slack nodes) connected. Parameters for the TR-slack are also given in Table 9, and it had the same on-load tap changer as TRF1 and TRF2.

The WPP park consists of 25 units of 5 MW, so the total generation was set to 125 MW and 0 Mvar.

The GFPP1 has an efficiency of 60%, and the operating point of active power was set to 10 MW, while the reactive power controlled the voltage of node B7 to 1.03 p.u.. The minimum allowed gas pressure was set to 85% of the nominal gas pressure at node EXIT1.

The GFPP2 also has an efficiency of 60%, and the minimum allowed gas pressure was 85% of nominal gas pressure at node EXIT6. The operating point was set to 70 MW and 2 Mvar.

Dissimilar to case study one, in case study two, the compressor driver consumption was set to 1.3% of the inlet gas mass flow to all compressors. For compressor CS2, which was the only electric-driven compressor, in this case, the driver consumption of 1.3% of the inlet gas mass flow was recalculated to the electric side via the above-described gas heating value at reference conditions. Pipe efficiency was not changed in this case, i.e., it was left as 100% (no additional gas pressure loss along the gas pipeline was introduced).

Within case study two, three scenarios (S7, S8, and S9) were simulated as described below. All simulations were performed for both UGS operation states: the injection and the withdrawal process.

The goal of scenario S7 was to identify the influence of the operation state of GFPP1, GFPP2, and P2G on the natural gas network node pressures and volumetric gas flows. The performed simulations showed that the same results were obtained for both the UGS operation processes, and they can be summarized as follows:

• When GFPP1 and GFPP2 were in operation, the most influenced gas nodes were N6 and N25, respectively, (gas pressures were decreased for 0.01 p.u. and 0.088 p.u., respectively). At the same time, the volumetric gas flows at nodes EXIT1 and EXIT6 were increased by 1.1% and 16.4%, respectively, both compared with situations when GFPP1 and GFPP2 were out of operation. The decreased gas pressures of N6 and N25 were still higher than the nominal gas pressures of nodes EXIT1 and EXIT6; therefore, the influences of GFPP1 and GFPP2 on the gas network could be neglected.
• When the P2G was in operation, the most influenced gas node was N12 (gas pressure was decreased by 0.014 p.u.). Interestingly, the volumetric gas flow produced by P2G (0.026 m³/s) was a little higher than the volumetric gas flow at node N14 demanded by EXIT3 (0.022 m³/s). This meant that the losses in the natural gas network could be reduced since natural gas did not need to flow anymore from node N9 to node N14 (130 km). Its flow reached velocities below the minimum allowed (less than 1.5 m/s). Thus, it can be concluded that P2G had a positive influence on the gas network.

In scenario S8, the influences of the CS1 operation state on the security of the power network supply were analyzed. In a normal operation, the CS1 regulates the output pressure to the nominal value. In the UGS injection process, when the CS1 was bypassed, the gas pressure of node N6 fell to 29.3 bar, as shown in Figure 5. Since the nominal pressure of node EXIT1 was 35 bar, GPRMS1 was bypassed, so gas pressures at nodes N6 and EXIT1 were equal. The gas pressure of 29.3 bar at node EXIT1 was lower than the minimum gas pressure required by GFPP1; thus, GFPP1 fell out of operation. This was reflected in the power network in such a way that lines B12-B13 and B9-B10 were overloaded (112.5% and 135.9%, respectively). The overloading of overhead lines in contingency situations is usually allowed up to 120%, so the overloading of line B12-B13 was not critical. The overloading of line B9-B10 was outside the upper limit; therefore, the load curtailment was expected in the 20 kV power network. This situation could be overcome if the UGS was in the withdrawal process. In that case, the gas pressure at node EXIT1 would be 0.89 p.u. and GPRMS1 would not fall out of operation.

In scenario S9, the influences of the CS2 operation state on the security of the power network supply, as well as the influences of a contingency in the power network on the CS2 driver power supply, were analyzed. The latter represents the influence of the power network on the security of the natural gas network supply. In a normal operation, CS2 regulates the output pressure to the nominal value. If CS2 was bypassed, the gas pressure at node N23 would become equal to the gas pressure at node N22, as shown in Figure 6. Consequently, in both UGS operation processes, the gas pressure at node EXIT6 would fall below the minimum required by GFPP2. Thus, GFPP2 would fall out of operation, which would cause voltage drops at B5, B6, B12, and B14 below the minimum allowed of 0.9 p.u., as well as the overloading of lines B1-B3, B3-B6, B9-B11, and B9-B10. Only the latter had an overload above the allowed 120%. Consequently, load curtailment in the power network could be expected to overcome this issue.

The same situation, as described in scenario S9, could occur if line B1-B2 fell out of operation. As a result, the CS2 driver would stay without a power supply, so the CS2 would have to be bypassed. Thus, the operation problems in the power network could be reflected in the natural gas network and returned to the power network, worsening the situation in the power network. Without the implementation of an integrated approach to modeling multi-energy systems, the problem mentioned above would not be observed. The proposed unique electrical model of a multi-energy system that uses the electrical analogy approach easily copes with this problem.

4. Conclusions

In this paper, the steady-state electrical equivalent models of the most important gas network elements (gas compressor station, gas pressure reducing and metering station, pipeline, gas load, gas source, and underground gas storage), as well as the steady-state electrical equivalent models of the linking facilities between the power and natural gas systems (gas-fired power plant and power-to-gas facility) were given. The equivalent models were developed using the network port theory and the load flow method formulation, known in the power systems analysis. In that way, the overall multi-energy system could be simultaneously solved as one electrical network and in one simulation step, i.e., without high resource-consuming iterative inter-model calculations associated with the separate network modeling approach. Moreover, the use of the extended Newton–Raphson method exceeded the calculation prerequisites associated with the ordinary Newton–Raphson method.

Two case studies were conducted. The developed models of network elements were used for modeling MES in the electrical software package—NEPLAN. In the first case study, the high accuracy of the presented models was confirmed by comparing the simulation results of one gas network modeled in the electrical analogy with the results obtained by SIMONE (a well-known software package for natural gas network simulations). The applicability of the presented approach was demonstrated in the second case study, where the electrical network was added to the gas network model together with the linking facilities. Finally, the case studies showed that the presented approach could be easily applied to the security of supply analyses (scenarios S8 and S9), as well as to determine the interdependencies between the coupled network infrastructures (scenario S7). Moreover, by conducting analyses on a multi-energy system, it was possible to detect some system vulnerabilities that could not be observed if the networks were modeled separately. For example, scenario S8 showed that operation problems in the electrical power network caused by outages in the natural gas network could be avoided if the operating regime of the underground gas storage was changed. Additionally, scenario S9 showed that the outages in the power network could affect the natural gas network, causing problems there, and, consequently, go back and worsen the situation in the power network.

The application of the presented methodology was limited to steady-state networks. Thus, to be able to use the model on an hourly basis, it is necessary to take into account natural gas stored in the gas pipelines, so-called linepack. The development of such a quasi-stationary model of a multi-energy system is the direction of our future research.