**Szymon Kuczy ´nski \*, Mariusz Łaciak, Andrzej Olijnyk, Adam Szurlej and Tomasz Włodek**

AGH University of Science and Technology, Drilling, Oil and Gas Faculty, 30-059 Krakow, Poland; laciak@agh.edu.pl (M.Ł.); aoliinyk@agh.edu.pl (A.O.); szua@agh.edu.pl (A.S.); twlodek@agh.edu.pl (T.W.) **\*** Correspondence: szymon.kuczynski@agh.edu.pl; Tel.: +48-12-617-2247

Received: 31 December 2018; Accepted: 20 February 2019; Published: 24 February 2019

**Abstract:** During the natural gas pipeline transportation process, gas stream pressure is reduced at natural gas regulation stations (GRS). Natural gas pressure reduction is accompanied by energy dissipation which results in irreversible exergy losses in the gas stream. Energy loss depends on the thermodynamic parameters of the natural gas stream on inlet and outlet gas pressure regulation and metering stations. Recovered energy can be used for electricity generation when the pressure regulator is replaced with an expander to drive electric energy generation. To ensure the correct operation of the system, the natural gas stream should be heated, on inlet to expander. This temperature should be higher than the gas stream during choking in the pressure regulator. The purpose of this research was to investigate GRS operational parameters which influence the efficiency of the gas expansion process and to determine selection criteria for a cost-effective application of turboexpanders at selected GRS, instead of pressure regulators. The main novelty presented in this paper shows investigation on discounted payback period (DPP) equation which depends on the annual average natural gas flow rate through the analyzed GRS, average annual level of gas expansion, average annual natural gas purchase price, average annual produced electrical energy sale price and CAPEX.

**Keywords:** natural gas; natural gas regulation station; turboexpander; pressure regulator; energy recovery; energy conversion; energy system analysis

#### **1. Introduction**

Natural gas is usually transported over long distances through pipelines at high pressures. In order to distribute the gas locally at points along the pipeline, the pressure must be significantly reduced before it is supplied to local distribution systems [1,2]. Currently, most pressure-reducing stations use expansion valves to reduce pressure [3,4]. When there is no heat transfer to or from the environment and if no work is done, the process of choked flow of the natural gas stream, regardless of its type, is an isenthalpic process [5]. In choked flow at constant enthalpy, the gas temperature changes. This change in gas temperature, once the pressure has been reduced, is associated with the so-called Joule–Thomson effect which has a negative value for methane which is the main component of natural gas (in the range of pressures and temperatures in the gas industry) [6]. For high-methane natural gas, it can be assumed with sufficient accuracy that the temperature drop equals approximately 0.4 ◦C/bar of pressure reduction [7]. While the gas flows through the station, the reduced natural gas temperature might produce hydrates inside the gas pipeline, ice plugs (especially in remotes for regulators) and icing of fittings, or contribute to the unstable functioning of a regulator work, etc. [8]. If the pressure regulator is replaced by expansion devices (expansion machine), the process will be isentropic (i.e., the system will perform external work, so that a significant part of the mechanical work will be recovered). However, in this case, the temperature drop is greater [9]. The choice of the type of expansion machine (gas expansion turbine or piston expansion device) which cooperates with a power generator depends

on the size of the GRS and gas flow parameters [10]. Due to this, it requires additional heating of natural gas which leads to additional costs. It is assumed that in order to generate 1 kWh of electricity by an expander, it is necessary to supply 1.15 kWh of heat [11]. The Polish natural gas transmission system has 1041 outlet points, and most of them are gas pressure regulation stations [12], but only a few dozen are suitable to use expansion machines and to produce economically viable electricity. The quantity of produced electricity depends on the gas flow rate through the gas regulation station (GRS) and the pressure reduction [13,14]. Some GRSs have large capacities, but the pressure reduction at these stations is quite low. On the other hand, there are small-capacity GRSs where the pressure can be reduced as much as 10 times relative to the inlet pressure. Such diversity in GRS technical parameters makes it difficult to select stations which could be equipped with expanders. Natural gas consumption in countries with high annual variability of the ambient temperature is highly irregular, which makes this selection even more difficult. This is because an expander's operational parameters deviate from nominal values, which has a negative impact on the efficiency of the device [12]. The economic impact which results from the application of expanders is also affected by the cost of the expander itself, as well as by the total investment cost [15]. Additionally, the sale price of produced electricity and the cost of natural gas heating has a large impact on the economic effect [16]. The purpose of this research was to develop an equation to assess the economic efficiency which results from the application of turboexpanders at individual stations.

Vortex turboexpanders usually have quite a high efficiency, in the range of 70–90%, which depends on the gas flow rate. Seasonal fluctuations in gas flow through the GRS cause an unequal generation of electricity by the turboexpander. Significant deviations in natural gas flow through the expander from the nominal numbers of flow indicated by the device manufacturer results in, among other things, a decrease in the efficiency of the device, and hence, less electricity production, even with large volumes of natural gas flow. In addition, a critically small or large natural gas flow through the expander can lead to its damage. Therefore, it is not recommended to use an expander as the only reducing device at the GRS. A typical diagram of the reduction and metering station is shown in Figure 1.

**Figure 1.** Scheme of the natural gas Reducing and Metering Station (without turboexpander). 1—Valve; 2—Filter; 3—Turbine gas meter; 4—Ultrasound gas meter; 5—Heat exchanger; 6—Shut down valve; 7—Pressure regulator; 8—Pressure drop valve; 9—Regulation valve; 10—Ventilation valve; 11—Gas odorizer.

A slightly more expensive, but much safer solution is to design a station which contains a traditional reduction sequence with an expander installed on it and a traditional reduction sequence which contains standard pressure reducers, as shown in Figure 2.

**Figure 2.** Schemes of Reducing and Metering Station (with turboexpander)**.** 1—Valve; 2—Filter; 3—Turbine gas meter; 4—Ultrasound gas meter; 5—Heat exchanger; 6—Shut down valve; 7—Pressure regulator; 8—Pressure drop valve; 9—Regulation valve; 10—Ventilation valve; 11—Gas odorizer; 12—Turboexpander; 13—Generator.

This solution ensures reliability and continuity of operation of the GRS as part of the natural gas transmission and distribution system. In situations where gas flow is too high or too low, the station's automatic control system or dispatching services will connect the reduction sequence, and the reduction will only be performed with a traditional reduction sequence. When natural gas flow stabilizes (back to the flow range recommended by the expander manufacturer), it will be possible to run the expander as a basic reduction element.

Due to the larger temperature drop caused by isentropic expansion, the reduction line which contains the expander should have a larger heat exchanger, or as shown in the diagram above, two heat exchangers in series. To ensure a stable pressure level after expansion, an ordinary regulator is recommended after the expander because gas pressure at the expander outlet depends on the gas pressure at the expander inlet. Application of the reducer can provide a stable pressure at the GRS outlet.

Issues related to the application of turboexpanders at GRS were discussed in scientific papers. Arabkoohsar et al., who presented a technical and economic investigation on the gas stations to check which type of power productive gas expansion stations are suitable to be accompanied with, combined heat and power systems [16]. The integration possibility of a pressure reduction station with low-temperature heat sources was studied by Borelli et al. [17]. They stated that the novel proposed system configuration has many advantages in terms of the opportunity to exploit low-enthalpy heat sources, highly efficient primary-conversion-technology utilizations, and integration with renewable

sources. Li et al. simulated the usage of a single-screw expander in the Jiyuanzhongyu natural gas letdown station in Henan province. The conclusion of their work was that the single-screw expander has great prospects in small-scale waste heat pressure energy recovery systems, as well as in new and renewable energy applications [18]. Stanek et al. presented the usefulness of exergy analysis, which in contrast to energy analysis, provides tools to indicate and quantify the origin of the electricity generated by a turboexpander [19]. Similar research was performed by Prilutskii, who compared the application of turbo- and piston expanders [20].

New technical solutions are also being developed to expand the possibilities of turboexpanders' application at small GRS. Barbarelli et al. developed a new microturbine typology which can operate with compressible fluids like steam or gases [21]. The novelty of this turbine occurs in its capability to operate with a low rotational speed with respect to the change in other operational parameters (i.e., flow rate, expansion ratio) without excessive loss of efficiency [22].

In most publications, the effectiveness of expander use is determined by economic indicators, whereas Lo Cascio et al. presented key performance indicators (KPIs) for integrated natural gas pressure reduction stations with energy recovery. Presented KPIs covered the different aspects of energy recovery in GRS. The waste energy recovery (WER) index is fundamental to assess the efficiency of the recovery process. Other indicators, such as the recovery ratio (RR) and the carbon emission recovery index (CER), provide information related to the status of the recovery process and carbon emission reduction, respectively [23]. In addition, GRS equipment optimization was discussed by Lo Cascio. Energy recovery from pressure reduction in natural gas networks was addressed through a structured retrofitting approach (SRA) aimed at optimal design. The SRA optimization technique permits the identification of the best system configuration in the GRS retrofit, particularly with regard to the turbo expander model and heat supplier technology. All the different issues of the integrated system design are addressed with the optimization technique [24].

#### **2. Calculation Model Description**

In order to develop selection criteria for gas regulation stations with turboexpanders, an equation was introduced including one factor which depends on five other independent parameters. This equation presents the impact of five independent parameters to determine the factor which indicates (or not) the application validity of the expander at the considered GRS. The discounted payback period (DPP) was chosen as a factor to present the economic efficiency of expanders in GRS. Some selected and discussed independent parameters are: annual average gas flow through the analyzed GRS, average annual level of gas expansion, average annual sale price of generated electrical energy, average purchase price of natural gas, and total investment cost. A mathematical model was developed to determine the impact of independent factors on the discounted payback period (DPP). Results of modeling were used as entry data for a statistical model (nonlinear polynomial regression). The calculation model consists of three sections: thermodynamic model, economic model, and statistical model.

#### *2.1. Thermodynamic Model*

Assuming that natural gas expansion in an expander is isentropic, the amount of electrical energy which can be produced can be calculated using the following equation [25]

$$\mathcal{W} = \eta\_{\mathcal{o}} \cdot \eta\_{\mathcal{c}\mathcal{S}} \cdot \eta\_{\mathcal{m}} \cdot M \cdot \Delta H \tag{1}$$

The gas enthalpy was calculated using [26,27]:

$$H = H\_0 + (G - G\_0) + T \cdot (S - S\_0) \tag{2}$$

The entropy of the gas was determined by using:

$$S - S\_0 = \frac{-\partial}{\partial T} (G - G^\sigma)\_p \tag{3}$$

Then, Gibbs energy was calculated by using:

$$\mathbf{G} - \mathbf{G}^o = \int\_{P\_0}^P V \cdot dP = \int\_0^P V \cdot dP + \int\_{P\_0}^0 V \cdot dP \tag{4}$$

where *V* = *<sup>R</sup>*·*<sup>T</sup> <sup>P</sup>*<sup>0</sup> , which after substitution will get,

$$\mathbf{G} - \mathbf{G}^o = \int\_0^P \left( V - \frac{\mathbf{R} \cdot \mathbf{T}}{P} \right) \cdot dP + \int\_{P\_0}^0 \left( V + \frac{\mathbf{R} \cdot \mathbf{T}}{P} \right) \cdot dP \tag{5}$$

$$G - G^{\circ} = \int\_{0}^{P} \left( V - \frac{R \cdot T}{P} \right) \cdot dP + R \cdot T \cdot \ln \frac{P}{P\_0} \tag{6}$$

By presenting the Clapeyron Equation as *Z* = *<sup>P</sup>*·*<sup>V</sup> <sup>R</sup>*·*<sup>T</sup>* <sup>→</sup> *<sup>V</sup>* <sup>=</sup> *<sup>Z</sup>*·*R*·*<sup>T</sup> <sup>P</sup>* we can write the following formula:

$$G - G^o = \int\_0^P \left(\frac{z \cdot R \cdot T}{P} - \frac{R \cdot T}{P}\right) \cdot dP + R \cdot T \cdot \ln\frac{P}{P\_0} \tag{7}$$

$$\mathbf{G} - \mathbf{G}^o = \mathbf{R} \cdot \mathbf{T} \int\_0^P \frac{\mathbf{1}}{P} (Z - \mathbf{1}) \cdot dP + \mathbf{R} \cdot \mathbf{T} \cdot \ln \frac{P}{P\_0} \tag{8}$$

$$G - G^{\circ} = R \cdot T \int\_{0}^{P} (Z - 1) \cdot d \ln P + R \cdot T \cdot \ln \frac{P}{P\_0} \tag{9}$$

Based on GERG 88 equation of state (EOS) [28],

$$Z = 1 + \frac{B}{R \cdot T} \cdot P + \frac{C - B^2}{\left(R \cdot T\right)^2} \cdot P^2 \tag{10}$$

where:

$$B = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \mathbf{x}\_{i} \cdot \mathbf{x}\_{j} \cdot B\_{ij}(T) \tag{11}$$

$$C = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \sum\_{k=1}^{n} \mathbf{x}\_{i} \cdot \mathbf{x}\_{j} \cdot \mathbf{x}\_{k} \cdot \mathbf{C}\_{ijk}(T) \tag{12}$$

$$B\_{ij}(T) = b\_{ij}^{(0)} + b\_{ij}^{(1)} \cdot T + b\_{ij}^{(2)} \cdot T^2 \tag{13}$$

$$\mathcal{L}\_{i\bar{j}\bar{k}}(T) = c\_{i\bar{j}k}^{(0)} + c\_{i\bar{j}k}^{(1)} \cdot T + c\_{i\bar{j}k}^{(2)} \cdot T^2 \tag{14}$$

By substituting Equation (10) into Equation (11) we obtain,

$$G - G^{\mathcal{O}} = R \cdot T \int\_{0}^{P} \left( \frac{B}{R \cdot T} \cdot P + \frac{C - B^2}{\left(R \cdot T\right)^2} \cdot P^2 \right) \cdot d \ln P + R \cdot T \cdot \ln \frac{P}{P\_0} \tag{15}$$

By solving Equation (15) we obtain the Gibbs energy:

$$G - G^o = RB \cdot P + \frac{\left(C - B^2\right) \cdot P^2}{2 \cdot R \cdot T} + R \cdot T \cdot \ln\frac{P}{P\_0} \tag{16}$$

From Equation (16), it is possible to derive the entropy equation. By substituting Equation (9) to Equation (3) we obtain

$$S - S\_0 = R \int\_0^P \left( 1 - Z - T \cdot \left( \frac{\partial Z}{\partial T} \right)\_P \right) \cdot d \ln P - R \cdot \ln \frac{P}{P\_0} \tag{17}$$

By substituting Equation (10) into Equation (17) we obtain:

$$S - S\_0 = R \int\_0^P \left( -\frac{B}{R \cdot T} \cdot P + \frac{C - B^2}{\left(R \cdot T\right)^2} \cdot P^2 - T \cdot \left(\frac{\partial Z}{\partial T}\right)\_P \right) \cdot d\ln P - R \cdot \ln \frac{P}{P\_0} \tag{18}$$

where:

$$
\left(\frac{\partial Z}{\partial T}\right)\_P = \frac{P}{R \cdot T} \cdot \left(\frac{dB}{bT} - \frac{B}{T}\right) + \frac{P^2}{\left(R \cdot T\right)^2} \cdot \left(\frac{d\mathbb{C}}{bT} - 2 \cdot B \cdot \frac{dB}{bT}\right) + \frac{2 \cdot P^2}{R^2 \cdot T^3} \cdot \left(B^2 - \mathbb{C}\right) \tag{19}
$$

Combining Equation (18) and (19) and integrating and simplifying the newly created dependence, we get,

$$S - S\_0 = \frac{-2 \cdot R \cdot T \cdot B \cdot P - 2 \cdot R \cdot T^2 \cdot P \cdot \left(\frac{dB}{BT} - \frac{B}{T}\right) - T \cdot P^2 \left(\frac{dC}{BT} - 2 \cdot B \cdot \frac{dB}{BT}\right) + P^2 \cdot \left(C - B^2\right)}{2 \cdot R \cdot T^2} - R \cdot \ln \frac{P}{P\_0} \tag{20}$$

By substituting Equations (16) and (20) into Equation (2) and then simplifying, we obtain the equation of natural gas enthalpy (21)

$$H = H\_0 + \frac{2 \cdot P^2 \cdot \left(\mathbb{C} - B^2\right) - 2 \cdot R \cdot T^2 \cdot P \cdot \left(\frac{dB}{BT} - \frac{B}{T}\right) - P^2 \left(\frac{dC}{BT} - 2 \cdot B \cdot \frac{dB}{BT}\right)}{2 \cdot R \cdot T} \tag{21}$$

where:

$$H\_0 = \int\_0^T \mathbb{C}\_{p0} \cdot dT \tag{22}$$

The coefficient of mechanical efficiency of the expander *η<sup>m</sup>* is calculated using Equation (23) [29,30]

$$\eta\_{\mathfrak{m}} = \left( A \cdot \left( \frac{\mathcal{Q}}{\mathcal{Q}\_{\mathfrak{n}}} \right)^2 + B \cdot \left( \frac{\mathcal{Q}}{\mathcal{Q}\_{\mathfrak{n}}} \right) + \mathcal{C} \right) \cdot \eta\_{\mathfrak{n}} \tag{23}$$

Equation (23) was obtained using an empirical method based on the average operating characteristic of turboexpanders. Figure 3 shows an example of the relationship between the efficiency *μ*/*μ<sup>n</sup>* of the expander and gas flow *Q*/*Qn*.

**Figure 3.** Relationship between the ratio of real and minimal efficiency vs. the ratio of real and nominal gas flow.

Apart from the generated electrical energy, the amount of heat used for heating up the gas before expansion is calculated using the thermodynamic model. The amount of heat needed for this purpose depends on the temperature drop due to the expansion. Gas expansion in an expander is an isotropic process, therefore, the drop in gas temperature will be much greater than in the case of normal pressure reduction and the associated Joule–Thomson effect. The gas temperature after expansion is determined using Equation (24):

$$T\_2 = T\_1 - \eta\_o \cdot (T\_1 - T\_{2ad}) \tag{24}$$

where *T*2*ad* is temperature of gas after expansion in the adiabatic process, calculated using Equation (25),

$$T\_{2ad} = T\_1 \cdot \left(\frac{P\_2}{P\_1}\right)^{\frac{k-1}{k}}\tag{25}$$

Due to technical and operational reasons, the gas temperature at the outlet of GRS cannot be lower than 273 K. Therefore, the gas temperature after heating is calculated using Equation (26),

$$T\_1 = \frac{273}{1 - \eta\_o \cdot \left(1 - \left(\frac{P\_2}{P\_1}\right)^{\frac{k-1}{k}}\right)}\tag{26}$$

The amount of natural gas needed to supply the necessary amount of heat (to preheat the gas before expansion) is calculated from Equation (27),

$$\mathcal{W}\_{hcut} = \frac{M \cdot \mathbb{C}\_p \cdot \Delta T}{3600 \cdot \eta\_{kgs\text{ boule}}} \tag{27}$$

where Δ*T* is the difference between gas temperature after heating, *T*1, and gas temperature at the inlet to GRS, *Tin*.

$$
\Delta T = T\_1 - T\_{\text{in}} \tag{28}
$$

In order to carry out a reliable economic analysis of the use of the expander, only the additional amount of heat needed to heat natural gas should be taken into account. The demand for additional heat is caused by the additional temperature drop compared to the use pressure reducer.

The temperature difference between the required gas temperature after heating, in the case of a turboexpander, and the gas temperature after heating, in the case of a pressure reducer, is given by the formula:

$$
\Delta T\_{add} = T\_1 - (4 \cdot \Delta P - T\_{in}) \tag{29}
$$

Hence, an additional amount of natural gas to supply additional heat, to preheat the gas is given by Equation (30),

$$\mathcal{W}\_{\text{heat\\_add}} = \frac{M \cdot \mathcal{C}\_p \cdot \Delta T\_{\text{add}}}{3600 \cdot \eta\_{gas\\_boiler}} \tag{30}$$

Based on the amount of produced electricity and amount of natural gas needed to provide the necessary amount of heat, the economic effect of using a turboexpander can be assessed.

#### *2.2. Economic Model*

The economic efficiency of the investment can be based on the discounted payback period (DPP) [31]:

$$DPP = \frac{CAPEX}{DCF\_t} \tag{31}$$

where discounted cash flow (DCF) is given by Equation (32) [32],

$$DCF\_{l} = \left[\frac{CF\_{1}}{\left(1+r\right)^{1}}\right] + \left[\frac{CF\_{2}}{\left(1+r\right)^{2}}\right] + \dots + \left[\frac{CF\_{l}}{\left(1+r\right)^{l}}\right] \tag{32}$$

Cash flow (CF) was calculated by using Equation (33) [33]:

$$CF = \left(\mathcal{W}\_{\text{ce}} \cdot EE\_{\text{price}} - \mathcal{W}\_{\text{heat}\_{\text{add}}} \cdot \mathcal{N}G\_{\text{price}}\right) - OPEX \tag{33}$$

The total capital expenditure includes the cost of design works and telemetry. Unexpected expenditures can be assessed as double the purchase cost of the device. The purchase cost of the turboexpander is determined using Equation (34) [34],

$$KZUI = \boldsymbol{\alpha} \cdot \boldsymbol{N}\_{el}^{-\beta} \tag{34}$$

Additional operational costs (OPEX) related to the operation and maintenance of the expander can be assumed to be 2% of the total investment cost. Another expenditure is the cost of natural gas used for heating up the gas before expansion. From the perspective of the operator, the difference in cash flow in the present state and after installation of the expander was also assumed. This means that the cost of heating up the gas in the existing system should be subtracted from the cost of heating it in the system equipped with an expander. These costs are paid by the operator, even though there is no income from the sales of electrical energy. In other words, only the cost of warming (above the level at which the Joule–Thomson compensation effect is reached in the regulators) was included.

The only positive element in the cash flow calculations will be the income from the sale of generated electrical energy.

#### *2.3. Statistical Model*

In order to determine the relationship between the discounted payback period and defined independent functions, a linear regression model was used, for which model coefficients were evaluated with the least squares method [35]. A function of dependent and independent parameters was described using Equation (35),

$$Y = X \cdot b + \varepsilon \tag{35}$$

In the least squares method, a vector of coefficients *b*, for which Equation (36) has minimum significance, is the solution of Equation (35).

$$F(b) = \sum\_{i=1}^{n} \varepsilon\_i^2 = \left(Y - X \cdot b\right)^T \cdot \left(Y - X \cdot b\right) = \varepsilon^T \cdot \varepsilon \tag{36}$$

Equation (37) is the necessity and sufficiency condition of minimum significance for function (36),

$$\frac{\partial F}{\partial b} = 2 \cdot X^T \cdot X \cdot b - 2 \cdot X^T \cdot Y = 0 \tag{37}$$

Hence Equation (38) is:

$$X^T \cdot X \cdot b = X^T \cdot Y \tag{38}$$

Matrix *XT*·*<sup>X</sup>* has the size *<sup>m</sup>* × *<sup>m</sup>* and the following structure:

$$\mathbf{X}^{T} \cdot \mathbf{X} = \begin{vmatrix} n & \sum \mathbf{x}\_{i1} & \cdots & \sum \mathbf{x}\_{ik} \\ \sum \mathbf{x}\_{i1} & \sum \mathbf{x}\_{i1}^{2} & \cdots & \sum \mathbf{x}\_{i1} \cdot \mathbf{x}\_{ik} \\ \vdots & \vdots & \vdots & \vdots \\ \sum \mathbf{x}\_{ik} & \sum \mathbf{x}\_{i1} \cdot \mathbf{x}\_{ik} & \cdots & \sum \mathbf{x}\_{ik}^{2} \end{vmatrix} \tag{39}$$

Vector *<sup>X</sup>T*·*<sup>Y</sup>* has m projections:

$$\mathbf{X}^T \cdot \mathbf{Y} = \begin{vmatrix} \sum y\_i \\ \sum y\_i \cdot x\_{i1} \\ \vdots \\ \sum y\_i \cdot x\_{ik} \end{vmatrix} \tag{40}$$

The solution obtained with the least squares method is defined with matrix expression (41):

$$b = \left(\mathbf{X}^T \cdot \mathbf{X}\right)^{-1} \cdot \left(\mathbf{X}^T \cdot \mathbf{Y}\right) \tag{41}$$

The efficiency of the regression was evaluated with a determination coefficient, *R*2. For a polynomial regression, the determination coefficient could be determined using Equation (42):

$$R^2 = 1 - \frac{\left(\boldsymbol{y} - \boldsymbol{X} \cdot \boldsymbol{b}\right)^T \cdot \left(\boldsymbol{Y} - \boldsymbol{X} \cdot \boldsymbol{b}\right)}{\left(\boldsymbol{Y} - \boldsymbol{Y}\right)^T \cdot \left(\boldsymbol{Y} - \boldsymbol{Y}\right)}\tag{42}$$

where - *y* is a vector of magnitude *n*, which consists of medium significance (43):

$$Y\_i = \frac{1}{n} \cdot \sum\_{i=1}^{n} Y\_i \tag{43}$$

#### **3. Entry Data for Calculation Model**

Real production data from fifteen operational GRS were used as entry data for the calculation model, as this should create a model which is closest to real-world conditions, and consequently, yield a higher accuracy of the searched dependence between independent factors and objective function.

Figures 4–7 present an hourly plot of gas flow intensity under normal conditions in the analyzed GRS as a function of time, as well as pressure level, at the inlet to the GRS for the selected four stations.

The plot covers a year of GRS operation. Summer minima and winter peaks of gas flow are associated with the seasonality of natural gas consumption in Poland and are clearly visible in the plot. The pressure at the GRS outlet is maintained at 0.3 MPa.

The economic analysis was based on the following assumptions:


Figures 4–7 show how GRS operating parameters differ during the year. Reduction and Metering Stations are mainly used in medium-pressure networks, which supply natural gas to municipal customers. The share of natural gas in the energy balance of electricity production in 2017 in Poland was 3.61% [37]. In autumn and winter, due to the lowering of the ambient temperature, natural gas consumption strongly increases, which is shown in Figures 4–7. However, when the ambient temperature rises in summer, the consumption of natural gas decreases significantly. Such strong fluctuations in GRS performance parameters mean that the process of turboexpander selection for individual GRS is complicated [38,39]. This does not allow for maintenance of a stable level of electricity production because by choosing an expander with a nominal capacity close to the maximum natural gas flow in autumn and winter, it will not be possible to produce enough electricity during the summer, since gas flow will be too low and this will reduce the efficiency of the expansion device [40,41].

**Figure 4.** Natural gas flow and pressure at regulation stations GRS #1 (before expansion) over a period of one year.

**Figure 5.** Natural gas flow and pressure at gas regulation stations (GRS) #2 (before expansion) over a period of one year.

**Figure 6.** Natural gas flow and pressure at regulation stations GRS #3 (before expansion) over a period of one year.

**Figure 7.** Natural gas flow and at pressure regulation stations GRS #5 (before expansion) over a period of one year.

The annual fluctuations in the GRS # 1 work parameters, presented in Figure 4, are the smallest among those presented. It is seen that in the summer months the daily irregularity of gas consumption during weekends decreases. The decrease in gas flow presented in Figure 4 is larger than presented in Figure 3 but not as drastic as in the case of GRS # 3 (operation profile is shown in Figure 5). In the range of 5000–6800 h, GRS # 3 almost doesn't function, and the expansion device does not produce electricity. The lack of gas flow through the station negatively affects the economic effect of the expander as the investment costs incurred for the purchase and installation of the expander device are not recovered.

The performance profile of GRS #5 presented in Figure 6, differs significantly from the GRSs discussed above. In this case, a significant increase in gas flow occurs in the range of 3800–5600 h, after which a significant reduction in gas flow is observed. The presented performance profile indicates that GRS #5 is one of the natural gas supply stations which provide gas to a ring-shaped gas distribution pipeline network in a large city. For the range of 3800–5600 h, all or almost all of the city's natural gas demand was covered by GRS #5. At the time of a dramatic drop in the flow (5600 h), the city's gas supply was most likely switched to other stations.

Another similar method to solve the problem of GRS flow decay during the summer is the ring system at the distribution network. A ring system is characteristic for large cities which are powered by several GRSs. In the winter, all of these GRSs work at maximum power, while in the summer, the gas flow to each of them is small. If it were possible to manage the gas flow in such a way that most of the city's gas demand would be covered by only one GRS, it would be an ideal place to use the turboexpander, since there would be no large seasonal fluctuations in gas flow through the station [42].

The pressure level for all presented stations within the analyzed range of time did not change drastically. Visible pressure fluctuations are caused by weekly unequal natural gas consumption. On non-working days, consumption is smaller, which means that the pressure in the network increases, whereas when consumption increases, the pressure in the network decreases. The ability to cover small

peaks of gas consumption is called the cumulative capacity of the transmission network. It results in a gas pressure increase in the transmission and distribution network during the night hours and a gas pressure reduction in the morning and evening hours. The transmission network also reacts with an increase in pressure on non-working days [43].

#### **4. Results of Calculations**

Based on data as described in the previous section, a computer simulation was performed to examine technical and economic effects related to the application of vortex turboexpanders at 15 selected Reduction and Metering Stations. Figures 8–15 present generated electricity, heat demand for natural gas heating, as well as expander operational efficiency calculated on the basis of Equation (22) for the selected (as described in the previous chapter) four GRSs.

For all stations, except GRS #5 in the range of 3800–5200 h, heat demand for natural gas heating before expansion is higher than the amount of generated electricity. This may seem unreasonable because more energy is lost than produced. However, this is because natural gas, which we burn for natural gas heating (before expansion), is a primary energy carrier and costs four times less than generated electricity [26]. For this reason, the application of expanders is reasonable from an economic point of view.

Figure 8 presents electricity generation and heat demand for GRS #1. The relationship between the amount of generated electricity and heat demand is linear. The increase in electricity generation, caused by the increase of gas flow through the expander, results in an increased heat demand for natural gas heating before expansion. The only deviation from this rule occurs in the time interval 1200–1300 h. In the indicated time range, heat demand increased despite the decline in electricity generation. Data analysis of GRS #1 operation profile, see Figure 4, showed that in the described time interval, gas flow increased simultaneously and inlet pressure was reduced. In the discussed time interval, there was a significant reduction in the efficiency of the turboexpander operation, see Figure 9.

Turboexpander efficiency reduction is caused by exceeding the nominal gas flow rate for which the efficiency of the expansion device is the largest.

A combination of several negative factors such as lowering of inlet pressure and exceeding nominal gas flow caused a reduction in electricity production. On the other hand, the increase in heat demand was due to the increase in gas flow, which needed to be heated before expansion [44].

A similar phenomenon was observed for GRS #2. The profile of electricity generation and heat demand is presented in Figure 10. Turboexpander efficiency changes are shown in Figure 11. For GRS #2, the phenomenon described above is more pronounced and occurs twice.

For the purpose of this research, the authors acquired real data which covered one year of operation for 15 natural gas reduction and measurement stations (GRS). The assessment of selected GRS economic efficiency depends on the determination of the discount payback period (DPP). DPP is calculated on the basis of cumulative annual electricity production and cumulative annual heat demand (heat carrier is natural gas consumed by GRS's own needs).

On the basis of real data which covered one year of operation for 15 natural gas reduction and measurement stations (GRS), it was possible to obtain 15 dependent functions *Yj* which describe the statistical model. Based on only 15 dependent functions *Yj* and 15 independent parameters *Xi* (*i* = 1, . . . ,5), it was not possible to determine the relationship between *Yj* and *Xi* with high accuracy.

**Figure 8.** Generated electrical energy and heat demand at GRS #1 over a period of one year.

7LPHRI\*56RSHUDWLRQKRXUV

**Figure 11.** Expander efficiency for GRS #2 over a period of one year.

**Figure 12.** Generated electrical energy and heat demand at GRS #3 over a period of one year.

**Figure 13.** Expander efficiency for GRS #3 over a period of one year.

**Figure 15.** Expander efficiency for GRS #5 over a period of one year.

In addition, the performance characteristics of GRS vary year to year, mainly due to, e.g., climate reasons (mild or severe winter) or change (increase or decrease) in the number of natural gas consumers.

Therefore, an attempt to estimate the dependence between selected parameters *Xi*, based on one year of measurements cannot ensure sufficient accuracy of the developed equation, *Yj*. Unfortunately, during the research, the authors did not get access to operational data of other GRSs. To increase the quantity of input data in the statistical model, data extension was applied to the operational parameters of each GRS.

Natural gas flow rate, GRS inlet pressure, and capital expenditures (CAPEX) data have been extended. Data extension of these parameters was performed with the creation of artificial values which were determined in steps of 10% in the range from −50% to +50% for each real value of a given parameter, *Xi*. Whereas, data of natural gas purchase price and generated electricity sell price have been extended in the range from −10% to +10% with a step of 2%.

As a result of data extension, an additional 750 samples were obtained. In total, 765 dependent functions *Yj* and 765 values of *Xi* parameters were used for analysis which was performed in the developed model. Samples for which DPP was greater than 40 years were removed. Finally, 607 data sets were used to develop the statistical analysis.

By application of the multiple linear regression method, formula (44) was obtained:

$$Y\_j = 45.227 - 10^{-4} \cdot X\_1 - 1.38 \cdot X\_2 + 0.0313 \cdot X\_3 - 0.0243 \cdot X\_4 + 3.114 \cdot 10^{-4} \cdot X\_5 \tag{44}$$

Formula (44) describes the relationship between discounted payback period (DPP) as a dependent parameter and selected independent parameters such as annual average gas flow rate (*X*1), average annual expansion level (*X*2), natural gas purchase price (*X*3), produced electrical energy sale price (*X*4), and capital expenditures (CAPEX) (*X*5). The coefficient of determination *R*<sup>2</sup> for the obtained equation (44) equals 0.543. The closer the coefficient of determination is to 1.00, the more accurate the regression model is. To improve the accuracy of the developed model, the influence of the individual independent parameter (*Xi*) on parameter *Yj* (DPP) was studied. Dependencies are presented in Figure 16.

Figure 16 (top) shows the dependence of the discounted payback period (DPP) vs. average annual gas flow through GRS *X*<sup>1</sup> (highest gas flow affects a faster return on investment) and DPP vs. average annual level of gas expansion *X*2. A higher level of gas expansion results in additional work which is performed by the expander, thus more electricity is generated which influences the faster reimbursement of costs incurred for the investment.

The relation presented in Figure 16 (middle) shows the dependence of DPP vs. purchase price of natural gas *X*<sup>3</sup> (the increase in the purchase price of natural gas adversely affects on DPP) and DPP vs. sale price of generated electrical energy *X*<sup>4</sup> (increase in the sale price of electricity produced positively affects the DPP).

The dependence of DPP vs. capital expenditures *X*<sup>5</sup> (CAPEX) is presented in Figure 15 (bottom). Usually low investment costs make the payback time of invested capital shorter. However, the presented results contradict this statement. Despite the fact that the price of a large expander is much higher than the price of a small one, a large turbo expander can operate at higher gas flow rates and can generate larger volumes of electricity. As a result, a properly selected and installed turboexpander on GRS can generate more electricity, which will be converted into positive cash flow.

Capital expenditures *X*<sup>5</sup> (CAPEX) depend on the installed electricity generation unit capacity. The amount of generated electricity is a function of the gas flow and expansion level. As a result, the dependence of DPP vs. capital expenditures *X*<sup>5</sup> (CAPEX) is similar to the dependency curve of DPP vs. average annual gas flow through GRS *X*<sup>1</sup> and DPP vs. average annual level of gas expansion *X*2.

None of the independent individual parameters *Xi* have a linear influence on the dependent parameter *Yj*, as shown in Equation (44).

&\$3(;86'

**Figure 16.** Discounted Payback Period (*Yj*) vs. individual independent parameter (*Xi*).

Because the influence of most factors on DPP is nonlinear, as shown in Figure 16, the statistical model should be modified. The previous form (44):

$$Y\_{\bar{\jmath}} = b\_0 + b\_1 \cdot X\_1 + b\_2 \cdot X\_2 + \dots + b\_5 \cdot X\_5 + \varepsilon \tag{45}$$

was replaced with (45):

$$Y\_{\hat{\jmath}} = b\_0 + b\_1 \cdot X\_1^{\hat{\varrho}\_1} + b\_2 \cdot e^{X\_2 \hat{\varrho}\_2} + \dots + b\_5 \cdot X\_5 + e \tag{46}$$

To restore the statistical model to the linear polynomial regression, the variable should be changed to: *<sup>Z</sup>*<sup>1</sup> <sup>=</sup> *<sup>X</sup>β*<sup>1</sup> <sup>1</sup> , *<sup>Z</sup>*<sup>2</sup> = *<sup>e</sup>X*2*β*<sup>2</sup> etc., after which a linear regression with new independent variables will be obtained [35]:

$$Y\_{\bar{j}} = b\_0 + b\_1 \cdot Z\_1 + b\_2 \cdot Z\_2 + \dots + b\_5 \cdot Z\_5 + e \tag{47}$$

By using the least-squares method, relations between the parameters *Xi* and *Yj* were determined, which are graphically shown in Figure 16. After the statistical analysis, Equation (48) was obtained, which is dependent on the discounted payback period (DPP) and annual average natural gas flow rate through the analyzed GRS (*X*1), average annual level of gas expansion (*X*2), average annual natural gas purchase price (*X*3), average annual produced electrical energy sale price (*X*4) and CAPEX (*X*5).

$$Y\_{\circ} = -52.91 + 1.1 \cdot 10^5 \cdot X\_1^{-0.965} + 65.34 \cdot X\_2^{-0.71} + 1.39 \cdot X\_3^{0.403} + 5.732 \cdot 10^5 \cdot X\_4^{-1.506} - 1.906 \cdot 10^9 \cdot X\_5^{-1.973} \tag{48}$$

The coefficient of determination *R*<sup>2</sup> for the obtained Equation (48) equals 0.768, which is much better than for the multiple linear model (44). The obtained Equation (48) allows for sufficient accuracy to validate the appropriateness of using turboexpanders on individual GRS. It is recommended to apply the above solution at the first stage of individual GRS selection for expansion installation. To make a final investment decision, a more detailed economic analysis, that was previously described with Equation (48), needs to be carried out for each selected station.

### **5. Discussion and Conclusions**

Due to high efficiency and relatively low cost, expansion machines can be widely used in natural gas transmission and distribution systems. Turboexpanders can be used where gas flow is present (average annual volumes of at least a few thousand normal cubic meters per hour) and the reduction of gas pressure takes place. This means that expansion machines can be used at control nodes, reduction and metering stations at industrial facilities, and reduction and metering stations at natural gas distribution networks, which are the most common elements of a natural gas network.

The main factor which negatively affects the economy of turboexpander application is seasonal unevenness in natural gas consumption (summer minima and winter peaks). When the gas flow significantly deviates from the nominal gas flow at which the expander's operational efficiency is the highest, the production of electricity drastically decreases due to the decrease in turboexpander efficiency. This negatively affects its economic efficiency, thus extending the payback period and reducing the cash flow. A solution to the above problem is possible if the distribution network operates in a ring system. In this case, in the summer, gas flow can be regulated and directed to GRS with an installed expander. In this way, it is possible to avoid gas flow reduction to the GRS with the expander. A lack of seasonal irregularities will positively affect the amount of electricity generation, and thus, will improve the economic efficiency of turboexpander application.

The choice of which GRS an expander should be installed on, and the assessment of its economic efficiency is complicated and requires knowledge of thermodynamics, mechanics, natural gas engineering, and economics.

The formula developed and presented in this article allows a quick and easy evaluation of the applicability of an expander on the selected GRS I◦. The discounted payback period (DPP) was the basis for the assessment of the investment profitability.

From a statistical point of view, the discounted payback period (DPP) is a dependent parameter, while independent parameters are the selected annual average gas flow rate (*X*1), average annual expansion level (*X*2), natural gas purchase price (*X*3), produced electrical energy sale price (*X*4), and capital expenditures (CAPEX) (*X*5).

The influence of independent parameters on the dependent parameter was analyzed. It was observed that this relationship was not linear. The equation developed by the multiple linear regression method had low accuracy, as determined by *R*2. The coefficient of determination *R*<sup>2</sup> for the obtained Equation (44) was 0.543.

After the implementation of non-linear multiple regression, Equation (48) was obtained, which is characterized by a much higher coefficient of determination *R*<sup>2</sup> = 0.768. The obtained Equation (48) allows sufficient accuracy to validate the appropriateness of using turboexpanders on a given individual GRS. It is recommended to use the above solution as the first step in the selection of the station for the installation of expansion machines. To make a final investment decision, it is necessary to perform a detailed economic analysis for each selected station which was previously examined with the developed Equation (48).

**Author Contributions:** Conceptualization, S.K., M.Ł. and A.O.; Formal analysis, S.K. and A.O.; Funding acquisition, M.Ł. and A.S.; Investigation, S.K. and A.O.; Methodology, A.O.; Supervision, M.Ł. and A.S.; Validation, S.K. and A.O.; Visualization, S.K., A.O. and T.W.; Writing—original draft, A.O.; Writing—review & editing, S.K. and T.W.

**Funding:** This work has received funding from the Statutory Research of Natural Gas Department at Drilling Oil & Gas Faculty, no. 11.11.190.555

**Conflicts of Interest:** The authors declare no conflict of interest.
