Computational Fluid Dynamics Methods for Gas Pipeline System Control

1. Introduction

At the present level of development of long, branched gas transmission networks (GTN), solving the problems of improving safety, efficiency and environmental soundness of operation of industrial pipeline systems calls for the application of methods of numerical simulation. The development of automated devices for technical inspection and process control, and availability of high-performance computer hardware have created a solid technical basis to introduce numerical simulation methods into the industrial practice of GTN analysis and operation. One of the promising approaches for numerical analysis of GTN operating is the development and application of high-accuracy computational fluid dynamics (CFD) simulators of modes of gas mixture transmission through long, branched pipeline systems (CFD-simulator) (Seleznev, 2007).

Actually, a CFD-simulator is a special-purpose software simulating, in “online” and “real time” modes with a high similarity and in sufficient detail, the physical processes of gas mixture transmission through a particular GTN. The development of a CFD-simulator focuses much attention to correctness of simulation of gas flows in the pipelines and to the impact produced by operation of relevant GTN gas pumping equipment (including gas compressor unit (GCU), valves, gas pressure reducers, etc.) and the environment upon the physical processes under study.

From the standpoint of mathematical physics, a CFD-simulator performs numerical simulation of steady and transient, non-isothermal processes of a gas mixture flow in long, branched, multi-line, multi-section gas pipeline network. Such simulation is aimed at obtaining high-accuracy estimates of the actual distribution (over time and space) of fluid dynamics parameters for the full range of modes of gas mixture transmission through the specific GTN in normal and emergency conditions of its operation, as well as of the actual (temporal) distribution of main parameters of GTN equipment operation, which can be expressed as functional dependencies on the specified controls on the GTN and corresponding boundary conditions. Theoretically, the high-accuracy of estimates of gas flow parameters is achieved here due to (Seleznev et al., 2005): (1) minimization of the number and depth of accepted simplifications and assumptions in the mathematical modeling of gas flows through long, branched, multi-section pipelines and gas compressor stations (CS) on the basis of adaptation of complete basic fluid dynamics models, (2) minimization of the number and depth of accepted simplifications and assumptions in the construction of a computational model of the simulated GTN, (3) improving methods for numerical analysis of the constructed mathematical models based upon results of theoretical investigation of their convergence and evaluation of possible errors of solution, (4) taking into account the mutual influence of GTN components in the simulation of its operation, (5) detailed analysis and mathematically formal description of the technologies and supervisor procedures for management of gas mixture transport at the simulated GTN, (6) automated mathematic filtration of occasional and systemic errors in input data, etc.

Input information required for work of a CFD-simulator is delivered from the Supervisory Control and Data Acquisition System (SCADA-system) operated at the simulated GTN.

CFD-simulator’s operating results are used for on-line control of the specific GTN, as well as in short-term and long-term forecasts of optimal and safe modes of gas mixture transport subject to fulfillment of contractual obligations. Also, a CFD-simulator is often used as base software for a hardware and software system for prevention or early detection of GTN failures.

For better illustration of the material presented in this chapter, but without loss of generality, further description of a CFD-simulator will be based on a sample pipeline network of a gas transmission enterprise. For the purpose of modeling, natural gas is deemed to be a homogenous gas mixture. A CFD-simulator of a gas transmission enterprise’s GTN is created by combining CS mathematical models into a single model of the enterprise’s pipeline system, by applying models of multi-line gas pipelines segments (GPS) (Seleznev et al., 2007). At that, in accordance with their process flow charts, the CS models are created by combining of GCU, dust catcher (DC) and air cooling device (ACD) models by applying mathematical models of connecting gas pipelines (CGP).

In a CFD-simulator, the control of simulated natural gas transmission through the GTN is provided by the following control commands: alteration of shaft rotation frequency of centrifugal superchargers (CFS) of GCU or their startup/shutdown; opening or closing of valves at a CS and valve platforms of multi-line GPS; alteration of the rates of gas consumption by industrial enterprises and public facilities; alteration of the gas reduction program at reduction units; alteration of the operation program at gas distributing stations; change in the program of ACD operating modes, etc. Therefore, simulated control in a CFD-simulator adequately reflects the actual control of natural gas transmission through pipeline networks of the gas transmission enterprise.

Generally, a CFD-simulator can be divided into three interrelated components (elements) (Seleznev et al., 2007). Each of these components is an integral part of the CFS-simulator.

The first system element is a computational scheme of a gas transmission enterprise pipeline system built on the basis of typical segments representing minimum distinctions from a comprehensive topology of an actual system considering the arrangement of valves, the system architecture, laying conditions, the process flow scheme of the system’s CS, etc. The second component is a database containing input and operative (current) data on time-dependent (owing to valves operation) system topology, pipeline parameters, process modes and natural gas transmission control principles for an actual gas transmission enterprise. The third component of a CFD-simulator is a mathematical software which operates the first two CFD-simulator elements.

The mathematical software includes (in addition to the computation core) a user interface environment imitating the operation of actual control panels located at gas transmission enterprises control centers in a visual form familiar to operators. This provides for faster training and, for the operator, easier adaptation to the CFD-simulator.

2. Simulation of multi-line GPS by CFD-simulator

Multi-line GPS are long, branched, multi-section pipelines. For numerical evaluation of parameters of steady and transient, non-isothermal processes of the gas mixture flow in multi-line GPS, a CFD-simulator uses a model developed by tailoring the full set of integral fluid dynamics equations to conditions of the gas flow through long branched pipeline systems. Transform of the 3D integral problem to an equivalent one-dimensional differential problem is implemented by accepting the minimum of required simplifications and projecting the initial system of equations onto the pipeline's geometrical axis. At that, a special attention is given to the adequacy of simulation of pipeline junction nodes where the 3D nature of the gas flow is strongly displayed.

In order to improve the accuracy of the simulations, it is reasonable to use the CFD-simulator in the “online” and “real time” modes for the numerical analysis of the given processes. There are two mathematical models of fluid flow through branched pipeline: heat-conductive model of pipeline junction and nonconductive model of pipeline junction. These models were developed by S. Pryalov and V. Seleznev at the turn of the century. These alternatives differ in a way of simulation of gas heat transfer within pipeline junction. The principle underlying the simulations is to observe the major conservation laws as strictly as possible. In practice the simultaneous implementation of the models makes it possible to find an accurate solution in short time.

The basis for the mathematical models of fluid flow through branched pipeline was the geometrical model of a junction (fig. 1) proposed by S. Pryalov (Seleznev et al., 2005). In this model, volume $${}^{(0)}V$$ can be depicted as a right prism with base area $$S_{base}$$ and height $H$ (see fig. 1а). For the prism lateral surface with linear dimensions $${}^{(n)}\delta$$ , true is the following relation: $${}^{(n)}\delta ={}^{(n)}f/H$$ , where $${}^{(n)}f$$ is the cross-sectional area of the pipe adjacent to the junction core $${}^{(0)}V$$ . It should be noted that the summarized volume of the joint is equal to $V={\displaystyle \underset{n=0}{\overset{N}{\cup }}{}^{\left(n\right)}V}$ , where $${}^{(n)}Vn=1\dots N$$ is the volume of an infinitely small section of the pipe adjacent to the junction core $${}^{(0)}V$$ (see fig. 1b). The prism base area can be represented as follows: $$S_{base}={\varsigma }_{base}{}^{(1)}\delta {}^2$$ , where $${\varsigma }_{base}$$ is the factor depending on the prism base geometry only. Now volume $${}^{(0)}V$$ can be determined by the following formula: $${}^{(0)}V=H{S}_{base}=$$ $$=H{\varsigma }_{base}{\left({}^{(1)}f/H\right)}^2={\varsigma }_{base}{}^{(1)}f{}^2/H$$ , which means that $$\underset{H\to \infty }{\lim }{}^{(0)}V=\underset{H\to \infty }{\lim }\left[{\varsigma }_{base}{}^{(1)}f{}^2/H\right]=0$$ . The application of this geometrical model enabled us to approximate compliance with mass, momentum and energy conservation laws at the pipelines junction.

Simplifications and assumptions used to construct the heat-conductive model of pipeline junction include the following: (1) when gas mixture flows join together, pressure can change with time, but at each time step it will have the same value at the boundaries of the pipeline junction, (2) the simulations take account of «downwind” heat and mass exchange due to heat conduction and diffusion, (3) in the pipeline junction, the gas mixture instantaneously becomes ideally uniform all over the pipeline junction volume $${}^{\left(0\right)}V$$ (see fig. 1b), (4) effects of gas mixture viscosity in the pipeline junction (inside the volume $${}^{\left(0\right)}V$$ ) can be ignored, (5) there are no heat sources in $${}^{\left(0\right)}V$$ (inside the volume $${}^{\left(0\right)}V$$ ), (6) pipeline diameters near the pipeline junction are constant.

Figure 1.

A schematic of a pipeline junction (а – 3D drawing; b – diagram)

Then, the heat-conductive fluid dynamics model of a transient, non-isothermal, turbulent flow of a viscous, chemically inert, compressible, heat-conductive, multi-component gas mixture through multi-line GPS which consist of pipes of round cross-sections and rigid rough heat-conductive walls is represented in the following way (Seleznev et al., 2005):

for each pipe adjacent to the junction node

$$\frac{\partial \left(\rho f\right)}{\partial t}+\frac{\partial }{\partial x}\left(\rho wf\right)=0$$ E1

$$\frac{\partial }{\partial t}\left(\rho Y_{m}f\right)+\frac{\partial }{\partial x}\left(\rho Y_{m}wf\right)-\frac{\partial }{\partial x}\left(\rho f{D}_m\frac{\partial Y_m}{\partial x}\right)=0m=\overline{1N_S-1}{Y}_{N_S}=1-{\displaystyle \sum _{m=1}^{N_S-1}{Y}_m}$$ E2

$$\frac{\partial \left(\rho wf\right)}{\partial t}+\frac{\partial \left(\rho w^{2}f\right)}{\partial x}=-f\left(\frac{\partial p}{\partial x}+g\rho \frac{\partial z_1}{\partial x}\right)-\frac{\pi }{4}\lambda \rho w\left|w\right|R$$ E3

$$\begin{array}{l}\frac{\partial }{\partial t}\left[\rho f(\epsilon +\frac{w^2}{2})\right]+\frac{\partial }{\partial x}\left[\rho wf(\epsilon +\frac{w^2}{2})\right]=-\frac{\partial }{\partial x}\left(pwf\right)-\rho wfg\frac{\partial z_1}{\partial x}-\\ -p\frac{\partial f}{\partial t}+Qf+\frac{\partial }{\partial x}\left(kf\frac{\partial T}{\partial x}\right)-\Phi \left(T{T}_{\text{am}}\right)+\frac{\partial }{\partial x}\left(\rho f{\displaystyle \sum _{m=1}^{N_S}{\epsilon }_m{D}_m\frac{\partial Y_m}{\partial x}}\right)\end{array}$$ E4

for each of the junction nodes

$$\frac{\partial \rho }{\partial t}+{\displaystyle \sum _{n=1}^N{}^{(n)}\Theta {}^{(n)}\left(\frac{\partial \left(\rho w\right)}{\partial x}\right)}=0$$ E5

$$\begin{array}{l}\frac{\partial \left(\rho Y_m\right)}{\partial t}+{\displaystyle \sum _{n=1}^N{}^{(n)}\left(\frac{\partial \left(\rho Y_{m}w\right)}{\partial x}\right){}^{(n)}\Theta }-{\displaystyle \sum _{n=1}^N{}^{(n)}\left[\frac{\partial }{\partial x}\left(\rho D_m\frac{\partial Y_m}{\partial x}\right)\right]{}^{(n)}\Theta }=0\\ m=\overline{1N_S-1}{}^{(n)}Y{}_{N_S}=1-{\displaystyle \sum _{m=1}^{N_S-1}{}^{(n)}Y{}_m}\end{array}$$ E6

$${}^{(n)}\left(\frac{\partial \left(\rho w\right)}{\partial t}\right)+{}^{(n)}\left(\frac{\partial \left(\rho w^2\right)}{\partial x}\right)=-{}^{(n)}\left(\frac{\partial p}{\partial x}\right)-g\rho {}^{(n)}\left(\frac{\partial z_1}{\partial x}\right)-0.25{}^{(n)}\left(\frac{\lambda \rho w\left|w\right|}{R}\right)n=\overline{1N}$$ E7

$$\begin{array}{l}\frac{\partial \left(\rho \epsilon \right)}{\partial t}+{\displaystyle \sum _{n=1}^N{}^{(n)}\left(\frac{\partial \left(\rho \epsilon w\right)}{\partial x}\Theta \right)}=-{\displaystyle \sum _{n=1}^N{}^{(n)}\left(p\frac{\partial w}{\partial x}\Theta \right)}+0.25{\displaystyle \sum _{n=1}^N{}^{(n)}\left(\frac{\lambda \rho {\left|w\right|}^3\Theta }{R}\right)}+\\ +Q+{\displaystyle \sum _{n=1}^N{}^{(n)}\left(\left[\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)\right]\Theta \right)}-{\displaystyle \sum _{n=1}^N{}^{(n)}\left(\frac{\Phi (T{T}_{\text{am}})}{f}\Theta \right)}+{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{m=1}^{N_S}{}^{(n)}\left(\left[\frac{\partial }{\partial x}\left(\rho {\epsilon }_m{D}_m\frac{\partial Y_m}{\partial x}\right)\right]\Theta \right)}}\end{array}$$ E8

$$\begin{array}{l}{}^{(n)}T={}^{(\xi )}T\epsilon ={}^{(n)}\epsilon ={}^{(\xi )}\epsilon {}^{(n)}\left({\epsilon }_m\right)={}^{(\xi )}\left({\epsilon }_m\right)\rho ={}^{(n)}\rho ={}^{(\xi )}\rho {}^{(n)}p={}^{(\xi )}p\\ {}^{(n)}k={}^{(\xi )}k{}^{(n)}\left(D_m\right)={}^{(\xi )}\left(D_m\right)Y_m={}^{(n)}\left(Y_m\right)={}^{(\xi )}\left(Y_m\right){}^{(n)}\left(z_1\right)={}^{(\xi )}\left(z_1\right)\\ \text{for any}n\xi \in \overline{1N}\text{and}m\in \overline{1N_S}\end{array}$$ E9

$${\displaystyle \sum _{n=1}^N{}^{(n)}\left(wfs\right)}=0{\displaystyle \sum _{n=1}^N{}^{(n)}\left(\frac{\partial T}{\partial x}fs\right)}=0{\displaystyle \sum _{n=1}^N{}^{(n)}\left(\frac{\partial Y_m}{\partial x}fs\right)}=0$$ E10

$${}^{(n)}s=\{\begin{array}{l}1if\left({}^{(0)}\overrightarrow{n}\cdot {}^{(n)}\overrightarrow{i}\right)0\\ -1if\left({}^{(0)}\overrightarrow{n}\cdot {}^{(n)}\overrightarrow{i}\right)0\end{array}{}^{(n)}\Theta =\frac{{}^{(n)}V}{V}=\frac{{}^{(n)}f{}_L{}^{(n)}\gamma }{{\displaystyle \sum _{k=1}^N\left({}^{(k)}f{}_L{}^{(n)}\gamma \right)}}0{}^{(n)}\Theta 1{\displaystyle \sum _{n=1}^N{}^{(n)}\Theta }=1$$ E11

equation of state (EOS) and additional correlations:

$$\begin{array}{l}p=p\left(\left\{S_{\text{mix}}\right\}\right)\epsilon =\epsilon \left(\left\{S_{\text{mix}}\right\}\right)k=k\left(\left\{S_{\text{mix}}\right\}\right){\epsilon }_m={\epsilon }_m\left(\left\{S_{\text{mix}}\right\}\right)D_m=D_m\left(\left\{S_{\text{mix}}\right\}\right)m=\overline{1N_S}\\ T_1=T_2=\dots =T_{N_S}=T\end{array}$$ E12

where $$\rho$$ is the density of the gas mixture; $$f$$ is the flow cross-sectional area of pipeline; $$t$$ is time (marching variable); $$x$$ is the spatial coordinate over the pipeline's geometrical axis (spatial variable); $$w$$ is the projection of the pipeline flow cross-section averaged vector of the mixture velocity on the pipeline's geometrical axis (on the assumption of the developed turbulence); $$Y_m$$ is a relative mass concentration of the m component of the gas mixture); $$D_m$$ is a binary diffusivity of component m in the residual mixture; $$N_S$$ is the number of components of the homogeneous gas mixture; $$p$$ is the pressure in the gas mixture; $$g$$ is a gravitational acceleration modulus; $$z_1$$ is the coordinate of the point on the pipeline's axis, measured, relative to an arbitrary horizontal plane, upright; $$\pi$$ is the Pythagorean number; $$\lambda$$ is the friction coefficient in the Darcy – Weisbach formula; $$R=\sqrt{f/\pi }$$ is the pipe's internal radius; $$\epsilon$$ is specific (per unit mass) internal energy of the gas mixture; $$Q$$ is specific (per unit volume) heat generation rate of sources; $$k$$ is thermal conductivity; $$T$$ is the temperature of gas mixture; $${\epsilon }_m$$ is specific (per unit mass) internal energy of the $$m$$ component; $$T_m$$ is the temperature of the $$m$$ component; $$N$$ is the number of pipes comprising one junction (see (5–11)); $${}^{(n)}s{}^{(n)}\Theta$$ are auxiliary functions (re $${}^{(0)}\overrightarrow{n}{}^{(n)}\overrightarrow{i}$$ see fig. 1b below); $$\left\{S_{\text{mix}}\right\}$$ is a set of parameters of gas mixture. Function $$\Phi \left(T{T}_{\text{am}}\right)$$ is defined by the law of heat transfer from the pipe to the environment and expresses the aggregate heat flow through the pipe walls along perimeter $$\chi$$ of the flow cross-section with area $$f$$ ( $$\Phi \left(T{T}_{\text{am}}\right)0$$ is cooling), $$T_{\text{am}}$$ is the ambient temperature. To denote the relationship of a value to the pipe numbered by $$n$$ , we use a parenthesized superscript on the left side of the value, e.g.: $${}^{(n)}\rho$$ . In equations (1–12), we use physical magnitudes averaged across the pipeline's flow cross-section. The set of equations (1–12) is supplemented by the boundary conditions and conjugation conditions. As conjugation conditions it is possible to specify boundary conditions simulating a complete rupture of the pipeline and/or its shutoff, operation of valves, etc.

As was stated above, the energy equations (4) and (8) comprise function $$\Phi \left(T{T}_{\text{am}}\right)$$ describing the heat exchange between the environment and natural gas in the course of its pipeline transmission. The space-time distribution of function $$\Phi \left(T{T}_{\text{am}}\right)$$ is defined, in the CFD-simulator, at specified time steps of the numerical analysis of parameters of the transient mode of gas transmission by solving a series of conjugate 2D or 3D problems of heat exchange between the gas flow core and the environment (Seleznev et al., 2007).

Simulation of steady processes of gas mixture flow through multi-line GPS is a less complicated task compared to (1–12). These models can be easily derived by simplifying the set of equations (1–12).

When simulating the gas mixture flow through real multi-line GPSs, one uses meshes with a spatial cell size of 10m to 10,000m. As a result, a smooth growth (or decrease) of temperature in a span of about $${10}^{-5}m$$ will be simulated as a temperature jump. Difference “upwind” equations make it possible to find a solution, which is quite adequate for the process at issue, almost without any impact on the convergence and accuracy of the resulting solution. On the contrary, the schemes that use the principles of central differences as applied to this process can yield solutions with difference oscillations. This may impair the accuracy of such simulations.

To overcome this drawback, S. Pryalov and V. Seleznev in 2008 suggested using the nonconductive model of pipeline junction. Downwind (and upwind) heat conduction and diffusion in this model are ignored. There will be a temperature discontinuity on gas mixture transition through the pipeline junction node, but this dependence will be monotone along each pipeline.

The list of additional simplifications in setting up the nonconductive model of pipeline junction includes only one item: the model ignores the downwind heat and mass exchange in the gas mixture due to heat conduction and diffusion. The temperature of the gas mixture at the outlet boundaries of inlet pipelines is defined only by the mixture flow parameters (mainly, by the temperature) inside these pipelines.

As there is no downwind heat transfer mechanism from the volume $${}^{\left(0\right)}V$$ to the inlet pipelines, the temperature of the gas mixture at these boundaries may differ from the temperature inside the volume $${}^{\left(0\right)}V$$ . On the other hand, the outlet boundaries of the inlet pipelines are also boundaries of the volume $${}^{\left(0\right)}V$$ . For this reason, it does not seem to be correct to say that the mixture is uniformly intermixed all over the volume.

Thus, the nonconductive model of pipeline junction of a transient, non-isothermal, turbulent flow of a viscous, chemically inert, compressible, multi-component gas mixture through multi-line GPS which consist of pipes of round cross-sections and rigid rough heat-conductive walls contains equations (1), (3), (7), (11) and the following equations:

for each pipe adjacent to the junction node

$$\frac{\partial }{\partial t}\left(\rho Y_{m}f\right)+\frac{\partial }{\partial x}\left(\rho Y_{m}wf\right)=0m=\overline{1N_S-1}{Y}_{N_S}=1-{\displaystyle \sum _{m=1}^{N_S-1}{Y}_m}$$ E13

$$\begin{array}{l}\frac{\partial }{\partial t}\left[\rho f(\epsilon +\frac{w^2}{2})\right]+\frac{\partial }{\partial x}\left[\rho wf(\epsilon +\frac{w^2}{2})\right]=\\ =-\frac{\partial }{\partial x}\left(pwf\right)-\rho wfg\frac{\partial z_1}{\partial x}-p\frac{\partial f}{\partial t}+Qf-\Phi \left(T{T}_{\text{am}}\right)\end{array}$$ E14

for each of the junction nodes

$$\begin{array}{l}{\displaystyle \sum _{n\in {ℂ}_{\text{IN}}}^{}\left\{\frac{\partial }{\partial t}\left[\rho \left(p_{\text{Joint}}{}^{\left(n\right)}T\left\{{}^{\left(n\right)}Y{}_m|m=\overline{1N_s}\right\}\right)\right]{}^{\left(n\right)}\Theta \right\}}+\\ +\frac{\partial }{\partial t}\left[\rho \left(p_{\text{Joint}}{T}_{\text{Joint}}\left\{Y_{m\text{Joint}}|m=\overline{1N_s}\right\}\right)\right]{\displaystyle \sum _{n\in {ℂ}_{\text{OUT}}}^{}{}^{\left(n\right)}\Theta }+\\ +{\displaystyle \sum _{n=1}^N\left\{\frac{\partial }{\partial x}\left[\rho \left(p_{\text{Joint}}{}^{\left(n\right)}T\left\{{}^{\left(n\right)}Y{}_m|m=\overline{1N_s}\right\}\right){}^{\left(n\right)}w\right]{}^{\left(n\right)}\Theta \right\}}=0\end{array}$$ E15

$${}^{\left(n\right)}\left(\frac{\partial }{\partial t}\left(\rho Y_{m}f\right)\right)+{}^{\left(n\right)}\left(\frac{\partial }{\partial x}\left(\rho Y_{m}wf\right)\right)=0m=\overline{1N_S-1}{Y}_{N_S}=1-{\displaystyle \sum _{m=1}^{N_S-1}{Y}_m}n\in {ℂ}_{\text{IN}}$$ E16

$$Y_{m\text{,}\text{Joint}}={\displaystyle \sum _{n\in {ℂ}_{\text{IN}}}{}^{\left(n\right)}\left(\rho \left|w\right|f{Y}_m\right)}/{\displaystyle \sum _{n\in {ℂ}_{\text{IN}}}{}^{\left(n\right)}\left(\rho \left|w\right|f\right)}m=\overline{1N_S}$$ E17

$$\begin{array}{l}{}^{(n)}\left(\frac{\partial }{\partial t}\left[\rho (\epsilon +\frac{w^2}{2})\right]\right)+{}^{(n)}\left(\frac{\partial }{\partial x}\left[\rho w(\epsilon +\frac{w^2}{2})\right]\right)=\\ =-{}^{(n)}\left(\frac{\partial \left(pw\right)}{\partial x}\right)-{}^{(n)}\left(\rho wg\frac{\partial z_1}{\partial x}\right)+{}^{(n)}Q-{}^{(n)}\left(\frac{\Phi \left({}^{(n)}T{}^{(n)}T{}_{\text{am}}\right)}{f}\right)n\in {ℂ}_{\text{IN}}\end{array}$$ E18

$$h_{\text{Joint}}={\displaystyle \sum _{n\in {ℂ}_{\text{IN}}}{}^{\left(n\right)}\left(\rho \left|w\right|hf\right)}/{\displaystyle \sum _{n\in {ℂ}_{\text{IN}}}{}^{\left(n\right)}\left(\rho \left|w\right|f\right)}n\in \{\begin{array}{l}{ℂ}_{\text{IN}}\in \overline{1N}if{}^{\left(n\right)}\left(ws\right)\ge 0\\ {ℂ}_{\text{OUT}}\in \overline{1N}if{}^{\left(n\right)}\left(ws\right)0\end{array}$$ E19

$$T_{\text{Joint}}=T\left(p_{\text{Joint}}{\epsilon }_{\text{Joint}}\left\{Y_{m\text{Joint}}|m=\overline{1N_s}\right\}\right){\rho }_{\text{Joint}}=\rho \left(p_{\text{Joint}}{T}_{\text{Joint}}\left\{Y_{m\text{Joint}}|m=\overline{1N_s}\right\}\right)$$ E20

$${\epsilon }_{\text{Joint}}=h_{\text{Joint}}-p_{\text{Joint}}/{\rho }_{\text{Joint}}{\displaystyle \sum _{n\in {ℂ}_{\text{IN}}}^{}{}^{\left(n\right)}\left(\rho \left|w\right|f\right)}={\displaystyle \sum _{n\in {ℂ}_{\text{OUT}}}^{}{}^{\left(n\right)}\left(\rho \left|w\right|f\right)}$$ E21

$$\begin{array}{l}{}^{(n)}p=p_{\text{Joint}}n=\overline{1N}{}^{(n)}\rho ={\rho }_{\text{Joint}}n\in {ℂ}_{\text{OUT}}{}^{(n)}T=T_{\text{Joint}}n\in {ℂ}_{\text{OUT}}{}^{(n)}\epsilon ={\epsilon }_{\text{Joint}}n\in {ℂ}_{\text{OUT}}\\ \end{array}$$ E22

$${}^{(n)}Y{}_m=Y_{m\text{,}\text{Joint}}n\in {ℂ}_{\text{OUT}}m=\overline{1N_S}{}^{(n)}(z_1)={}^{(\xi )}(z_1)\text{for any}n\xi \in \overline{1N}$$ E23

equation of state:

$$h=\epsilon +p/\rho p=p\left(\rho T\left\{Y_m|m=\overline{1N_s}\right\}\right)\rho =\rho \left(pT\left\{Y_m|m=\overline{1N_s}\right\}\right)$$ E24

$$\epsilon =\epsilon \left(pT\left\{Y_m|m=\overline{1N_s}\right\}\right)T=T\left(p\epsilon \left\{Y_m|m=\overline{1N_s}\right\}\right)$$ E25

where the subscript “ $\text{Joint}$ ” means that the physical parameter of the gas mixture belongs to the pipeline junction i.e. to the volume $${}^{\left(0\right)}V$$ $$h_{\text{Joint}}$$ is enthalpy of the gas mixture in the pipeline junction.

The numerical analysis of the mathematical models (1–12) and (1, 3, 7, 11, 13–24) under consideration will be carried out by hybrid algorithm. It comprises Integro-Interpolation Method by A. Tikhonov and A. Samarsky (IIM) (Russian analog of the Finite Volume Method) (Tikhonov & Samarsky, 1999) and Lagrangian Particle Method by Pryalov (LPM).

To illustrate the parametric classes used for the difference equations in IIM, it is possible to present the class of the difference equations for a mathematical model of the non-isothermal transient motion of a multi-component gas mixture through a GPS line (see (1–4, 12)) (Seleznev, 2007):

$${\left[{\left(\rho F\right)}^{\left(-S\right)}\right]}_{\underset{+}{t}}+{\left({\rho }_{\left(-R\right)}{w}_{\left(0.5\right)}{f}_{\left(0.5\right)}\right)}_{\underset{+}{x}}^{\left(\sigma \theta \right)}=0$$ E26

$$\begin{array}{l}{}^{}\left[{\left(\rho F\right)}^{(-S)}{\left(Y_m\right)}^{(-S)}\right]{}_{\underset{+}{t}}+\left({\rho }_{(-R)}{w}_{(0.5)}{f}_{(0.5)}{\left(Y_m\right)}_{(-R)}\right){}_{\underset{+}{x}}^{(\sigma \theta )}-\\ -\left[{\left(\rho f{D}_m\right)}_{(-R)}\delta {\left(Y_m\right)}^a\right]{}_{\underset{+}{x}}^{(\sigma \theta )}=0m=\overline{1N_S-1}{Y}_{N_S}=1-{\displaystyle \sum _{m=1}^{N_S-1}{Y}_m}\end{array}$$ E27

$$\begin{array}{l}{}^{}\left[{\left(\rho F\right)}^{(-S)}{w}^{(-S)}\right]{}_{\underset{+}{t}}+\left({\rho }_{(-R)}{w}_{(0.5)}{f}_{(0.5)}{w}_{(-R)}\right){}_{\underset{+}{x}}^{(\sigma \theta )}=\\ =-\left(B^{-}{\gamma }^{-}{p}_{\overline{x}}+B^{+}{\gamma }^{+}{p}_x\right){}^{(\sigma \theta )}-g\left[\rho \left(B^{-}{\gamma }^{-}{\left(z_1\right)}_{\overline{x}}+B^{+}{\gamma }^{+}{\left(z_1\right)}_x\right)\right]{}^{(\sigma \theta )}-\frac{\pi }{4}\left(\lambda \rho \left|w\right|r\right){}^{(\sigma \theta )}{w}^{(\sigma \theta )}\end{array}$$ E28

$$\begin{array}{l}\left[{\left(\rho F\right)}^{(-S)}{\epsilon }^{(-S)}\right]{}_{\underset{+}{t}}+\left({\rho }_{(-R)}{w}_{(0.5)}{f}_{(0.5)}{\epsilon }_{(-R)}\right){}_{\underset{+}{x}}^{(\sigma \theta )}+K_t\left(\rho F\frac{w^2}{2}\right)+K_x\left(\rho wf\frac{w^2}{2}\right){}^{(\sigma \theta )}=\\ =-\left(p_{(-R)}{w}_{(0.5)}{f}_{(0.5)}\right){}_{\underset{+}{x}}^{(\sigma \theta )}-g\left[\rho \left(B^{-}{\gamma }^{-}{\left(z_1\right)}_{\overline{x}}+B^{+}{\gamma }^{+}{\left(z_1\right)}_x\right)\right]{}^{(\sigma \theta )}{w}^{(\sigma \theta )}-p^{(\sigma \theta )}\left[F^{(-S)}\right]{}_{\underset{+}{t}}+\\ +{\left(QF\right)}^{(\sigma \theta )}+\left[{\left(kf\right)}_{(-R)}\delta T^a\right]{}_{\underset{+}{x}}^{(\sigma \theta )}-{\varphi }^{(\sigma \theta )}+{\displaystyle \sum _{m=1}^{N_S}\left[{\left(\rho D_{m}f\right)}_{(-R)}{\left({\epsilon }_m\right)}_{(-R)}\delta {\left(Y_m\right)}^a\right]{}_{\underset{+}{x}}^{(\sigma \theta )}}\end{array}$$ E29

$${\epsilon }_m={\epsilon }_m\left(\left\{S_{\text{mix}}\right\}\right)m=\overline{1N_S}{T}_1=T_2=\dots =T_{N_S}=Tp=p\left(\left\{S_{\text{mix}}\right\}\right)\epsilon =\epsilon \left(\left\{S_{\text{mix}}\right\}\right)$$ E30

$$k=k\left(\left\{S_{\text{mix}}\right\}\right)D_m=D_m\left(\left\{S_{\text{mix}}\right\}\right)m=\overline{1N_S}$$ E31

where $$F$$ $$B^{+}$$ $$B^{-}$$ are the expressions approximating $$f$$ $$r$$ is the expression approximating $$R$$ (the type of these expressions is defined upon selection of a particular scheme from the class of schemes); $$s^a{s}^b{r}^a{r}^b$$ $$\sigma$$ $$\theta$$ are parameters of the class of schemes (e.g., by specifying $$s^a=s^b=1r^a=r^b=05$$ $$\sigma =1$$ $$\theta =0$$ , a two-layer scheme with central differences is selected from the class of scheme, and by specifying $$s^a=s^b=1r^a{r}^b$$ is according to the principles of “upwind” differencing, $$\sigma =1$$ $$\theta =0$$ is a two-layer “upwind” scheme); $$K_t$$ and $$K_x$$ are the differential-difference operators of functions $$\left(0.5\rho F{w}^2\right)$$ and $$\left(0.5\rho wf{w}^2\right)$$ over time and space, respectively (the type of these operators is defined upon selection of a particular scheme from the class of schemes); $${\varphi }^{(\sigma \theta )}$$ is a difference expression approximating function $$\Phi \left(T{T}_{\text{am}}\right)$$ . The difference equations (25–30) are supplemented by difference expression of initial and boundary conditions, as well as of conjugation conditions.

To record the parametric class of the difference equations (25–30), we used notations of a non-uniform space-time mesh $${\left\{x_i{t}_j\right\}}^{}$$ , where $$x_i$$ and $$t_j$$ are coordinates of the mesh node numbered $$i$$ over space, and $$j$$ , over time, $$ij\in Z$$ $$Z$$ being a set of nonnegative integers.

To explain the notations, it is expedient to consider an individual computational cell containing the node $$\left\{x_i{t}_j\right\}$$ (mesh base node) and bounded by straight lines $$x=x_i^a$$ $$x=x_i^b$$ $$t=t_j^a$$ and $$t=t_j^b$$ $$x_{i-1}\le x_i^a\le x_i$$ $$x_i\le x_i^b\le x_{i+1}$$ $$t_{j-1}\le t_j^a\le t_j$$ $$t_j\le t_j^b\le t_{j+1}$$ $$x_i^a\ne x_i^b$$ $$t_j^a\ne t_j^b$$ Let us introduce the so-called weighing parameters:

$$r_i^a=\frac{x_i^a-x_{i-1}}{x_i-x_{i-1}}=\frac{x_i^a-x_{i-1}}{h_i}{r}_i^b=\frac{x_i^b-x_i}{x_{i+1}-x_i}=\frac{x_i^b-x_i}{h_{i+1}}{s}_j^a=\frac{t_j^a-t_{j-1}}{t_j-t_{j-1}}=\frac{t_j^a-t_{j-1}}{{\tau }_j}{s}_j^b=\frac{t_j^b-t_j}{{\tau }_{j+1}}$$ E32

where $$h_i=x_i-x_{i-1}$$ and $$h_{i+1}=x_{i+1}-x_i$$ are steps “backward” and “forward” over the space coordinate for the i node; $${\tau }_j=t_j-t_{j-1}$$ and $${\tau }_{j+1}=t_{j+1}-t_j$$ are steps “backward” and “forward” over the time coordinate on the $$j$$ time layer; $${\alpha }_i=h_{i+1}/h_i$$ and $${\beta }_j={\tau }_{j+1}/{\tau }_j$$ are mesh parameters characterizing non-uniformity of the space and time mesh. To describe mesh function $$y=y(xt)$$ , the system (25–30) used the following notations:

$$y_i^j=y\left(x_i{t}_j\right)y_i^{}=y_i^{}\left(t\right)=y\left(x_{i}t\right)y_{}^j=y_{}^j\left(x\right)=y\left(x{t}_j\right)$$ E33

Where there was applied the quadratic approximation

$$\begin{array}{l}y\left(xt\right)=a_i^{y-}\left(t\right)x^2+b_i^{y-}\left(t\right)x+c_i^{y-}\left(t\right)x\in \left[x_{i-1}{x}_i\right]\\ y\left(xt\right)=a_i^{y+}\left(t\right)x^2+b_i^{y+}\left(t\right)x+c_i^{y+}\left(t\right)x\in \left[x_i{x}_{i+1}\right]\end{array}$$ E34

the system used the following notations:

$$\begin{array}{l}\delta y\left(xt\right)=\frac{\partial y}{\partial x}=2a_i^{y-}(t)x+b_i^{y-}(t)x\in \left[x_{i-1}{x}_i\right]\delta y\left(xt\right)=\frac{\partial y}{\partial x}=2a_i^{y+}(t)x+b_i^{y+}(t)x\in \left[x_i{x}_{i+1}\right]\\ {\left(\delta y_{}^a\right)}_i^j=2a_i^{y-}\left(t_j\right)x^a+b_i^{y-}\left(t_j\right){\left(\delta y_{}^b\right)}_i^j=2a_i^{y+}\left(t_j\right)x^b+b_i^{y+}\left(t_j\right)\end{array}$$ E35

Also, the system (25–30) used the following index-free notations:

$$\begin{array}{l}y=y_i^{j}h=h_i\alpha ={\alpha }_i\tau ={\tau }_j\beta ={\beta }_{j}x=x_{i}t=t_j\stackrel{⌢}{y}=y_i^{j+1}\stackrel{⌣}{y}=y_i^{j-1}y(+1)=y_{i+1}^j\\ y(-1)=y_{i-1}^j{r}^a=r_i^a{r}^b=r_i^b{s}^a=s_j^a{s}^b=s_j^{b}ij\in Z{y}^{\left(a\right)}=ay+(1-a)\stackrel{⌣}{y}\\ y_{(b)}=by+(1-b)y(-1)y^{(ab)}=a\stackrel{⌢}{y}+(1-a-b)y+b\stackrel{⌣}{y}{y}^{(-S)}=s^{a}y+(1-s^a)\stackrel{⌣}{y}\\ y^{(+S)}=s^b\stackrel{⌢}{y}+(1-s^b)y{y}_{\left(-R\right)}=r^{a}y+\left(1-r^a\right)y\left(-1\right)y_{\left(+R\right)}=r^b\cdot y\left(+1\right)+\left(1-r^b\right)y\\ y_{\left(-R\right)}\left(+1\right)=y_{\left(+R\right)}{y}_{\left(+R\right)}\left(-1\right)=y_{\left(-R\right)}{\stackrel{⌢}{y}}^{\left(-S\right)}=y^{\left(+S\right)}{\stackrel{⌣}{y}}^{\left(+S\right)}=y^{\left(-S\right)}\Delta t^{-}=t-t^a=(1-s^a)\tau \\ \Delta t^{+}=t^b-t=s^b\stackrel{⌢}{\tau }=\beta s^b\tau \Delta x^{-}=x-x^a=(1-r^a)h\Delta x^{+}=x^b-x=r^{b}h(+1)=\alpha r^{b}h\\ \Delta t=t^b-t^a=\Delta t^{+}+\Delta t^{-}=(1-s^a)\tau +s^b\stackrel{⌢}{\tau }=(1-s^a)\tau +\beta s^b\tau \\ \Delta x=x^b-x^a=\Delta x^{+}+\Delta x^{-}=(1-r^a)h+r^{b}h(+1)=(1-r^a)h+\alpha r^{b}h\\ {\gamma }^{-}=\Delta x^{-}/\Delta x=(1-r^a)h/\left[(1-r^a)h+r^{b}h(+1)\right]=(1-r^a)/\left[(1-r^a)+\alpha \cdot r^b\right]\\ {\gamma }^{+}=\Delta x^{+}/\Delta x=r^{b}h(+1)/\left[(1-r^a)h+r^{b}h(+1)\right]=\alpha \cdot r^b/\left[(1-r^a)+\alpha \cdot r^b\right]\\ \delta y^a(+1)=\delta y^b\delta y^b(-1)=\delta y^a{y}_x=\left[y(+1)-y\right]/h(+1)y_{\overline{x}}=\left[y-y(-1)\right]/h\\ y_{\overline{x}}(+1)=y_x{y}_x(-1)=y_{\overline{x}}{y}_{\underset{+}{t}}=\left[\stackrel{⌢}{y}-y\right]/\Delta t{y}_{\underset{+}{x}}=\left[y(+1)-y\right]/\Delta x\end{array}$$ E36

LPM is illustrated by the example of solution of energy equation from gas dynamics equations set for a single-component gas transmitted through an unbranched pipeline. Imaginary Lagrange particles are distributed along the pipeline. They are considered weightless. This allows them to move together with the fluid. Due to the small size, each particle can instantaneously acquire the temperature of the ambient fluid. Thus, by tracking the motion of such Lagrange particles along with the fluid and their temperature, one can analyze the process of heat transfer through multi-line GPSs. Energy equation is easy to derive by simplifying and transforming equation (4) accounting for (1-3), (12) and (23). The simplified equation will have the following form:

$$\frac{\partial h}{\partial t}+w\frac{\partial h}{\partial x}=\frac{w}{\rho }\frac{\partial p}{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial t}+\frac{\lambda {\left|w\right|}^3}{4R}-\frac{\Phi \left(T{T}_{\text{am}}\right)}{f\rho }\text{or} \frac{Dh}{Dt}=G\left(xtT\left(h\right)\right)$$ E37

where

$$G\left(xtT\right)=\frac{w}{\rho }\frac{\partial p}{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial t}+\frac{\lambda {\left|w\right|}^3}{4R}-\frac{\Phi \left(T{T}_{\text{am}}\right)}{f\rho }$$ E38

$$D\xi /Dt=\partial \xi /\partial t+w\partial \xi /\partial x$$ is a derivative of an arbitrary function $\xi$ over $t$ in the direction

$$dx/dt=w\left(xt\right)$$ E39

This direction is called characteristic, and the equation is called the equation of characteristic direction. The second equation in (31) is called the characteristic form of the first equation in (31) or the differential characteristic relation. From the physical standpoint, the derivative $$Dh/Dt$$ corresponds to the substantial derivative, and the solution of equation (33) defines the coordinate of the continuum particle (in our case, the spatial coordinate of the fluid flow cross section) at each time.

Considering the known thermodynamic relationship

$$dh=c_{p}dT-\mu c_{p}dp$$ E40

where $c_p$ is specific heat capacity at constant pressure; $\mu$ is the Joule-Thomson factor, equation (31) can be transformed as follows:

$$\frac{DT}{Dt}=\left(\mu +\frac{1}{\rho c_p}\right)\frac{Dp}{Dt}+\frac{\lambda {\left|w\right|}^3}{4R{c}_p}-\frac{\Phi \left(T{T}_{\text{am}}\right)}{f\rho c_p}$$ E41

Equation (35) is satisfied along each characteristic curve (33). These curves per se describe the trajectory of fluid particles. In other words, these equations describe the change in the fluid temperature for each cross section of the transported product flow.

When implementing the LPM, fluid flow parameters (such as pressure and velocity) are obtained using a difference scheme, while the gas temperature distribution is obtained based on the analysis of the Lagrange particle motion. For each particle, equation (35) is solved. The form of this equation enables such simulations, because it corresponds to the change in time of the temperature at each cross section of the fluid flow. Within this problem statement (which is Lagrangian with respect to each particle), equation (35) transforms into an ordinary differential equation (ODE) relative to the marching variable:

$$\frac{dT}{dt}-\left(\mu +\frac{1}{\rho c_p}\right)\frac{dp}{dt}=\frac{\lambda {\left|w\right|}^3}{4R{c}_p}-\frac{\Phi \left(T{T}_{\text{am}}\right)}{f\rho c_p}$$ E42

Numerical analysis of ODE (36) can be carried out using different ODE solution procedures, for example, the known Runge – Kutta method with an adjustable accuracy of solution. As initial temperature of Lagrange particles one uses respective values from the defined initial conditions (i.e., for each Lagrange particle, its temperature is assumed equal to the fluid temperature at the location of the given particle).

As particles move together with the fluid flow towards the outlet pipe boundary, one needs to introduce new Lagrange particles at the inlet boundary at some regular intervals. The initial temperature of each introduced particle should be defined based on the boundary conditions related to the inlet boundary of the given pipe. The Lagrange particles that leave the pipe are deleted. As applied to the inlet boundaries of the outlet pipelines of each junction node, the temperature of the introduced particles should be defined in accordance with equations (10) and (12).

Since LPM for solving the energy equation is not related immediately to the finite difference mesh used to solve the continuity and momentum equations, this mesh has almost no effect on the accuracy of the method proposed. Thus, high-accuracy computed values of gas temperature are obtained without mesh refinement, and this leads to a considerable speedup of computations.

Also, as there is no direct relationship between LMP and the finite difference mesh, the method is free of artificial viscosity and computational dispersion (Fletcher, 1988). As a result, LPM yields solutions without “scheme smoothing” of temperature fronts, which is consistent with actual physical processes. This makes such simulations more credible than the simulations, in which the energy equation is solved using difference schemes.

3. Simulation of a CS by CFD-simulator

The principal task of mathematical simulation of stable and safe operation of a CFS is to determine physical parameters of gas at the CFS outlet on the basis of the known values of gas flow parameters at the CFS inlet. To construct a 1D mathematical model of a CFS in a CFD-simulator, we used a well-known polytropic model of a CFS developed by A. Stepanov. The model is based on the combination of analytical dependencies for polytropic fluid dynamics processes and empirical characteristics obtained for each CFS during its full-scale testing.

When simulating steady modes of CS operation, an isothermal model is used for description of the gas flow in a CGP and DC, and an isobaric model – for description of the gas flow in an ACD. The power drive is simulated by specifying an analytical dependency of the capacity at the CFS shaft on energy expenditures. Such approach provides for the simplicity of the conjugation of models and a high, from the practical standpoint, veracity of simulation.

As was noted above, a CFD-simulator of a particular CS is a result of combining GCU, ACD and DC models, by application of CGP models, into a single integral network model of the CS in accordance with the process flow charts of the actual CS (fig. 2).

As proposed by V. Seleznev (Kiselev et al., 2005), in order to determine parameters of steady modes of natural gas transmission through a CS, generally, it is necessary to solve a system of nonlinear algebraic equalities under simple constraints on unknown variables. The system includes the law of conservation of gas masses at CGP branching points and one of the group of equations representing either conditions for conservation of mass flow rate at inlet and outlet CGPs in one branch, or conditions for equality of natural gas heads in parallel branches, where a branch is a segment of a pipeline system, which comprises an inflow (inlet) CGP, CFS and an outflow (outlet) CGP (see fig. 2).

As independent decision variables we used fractions of a mass flow rate of natural gas transmitted through separate branches of a CS, ratios of compression by compressor shops and ratios of compression by GCUs working as the first stage of transmitted gas compression at compressor shops. Such a set of variables allowed to reduce the problem dimension and narrow the range of search for problem solution, by more accurate specification of constraints on variables. This allowed to considerably save running time.

Figure 2.

The scheme of decision independent variables assignment for on-line technological analysis of gas transmission through CS

The mathematical model for the CS scheme presented in fig. 2 can be written as follows:

$$\{\begin{array}{l}{P}_{11}^1\left(J_{in}{P}_{in}{T}_{in}{J}_1\left(X_1\right)J_{11}\left(X_2\right)X_6\right)-P_{12}^1\left(J_{in}{P}_{in}{T}_{in}{J}_1\left(X_1\right)J_{12}\left(X_3\right)X_4\right)=0\\ P_{12}^1\left(J_{in}{P}_{in}{T}_{in}{J}_1\left(X_1\right)J_{12}\left(X_3\right)X_4\right)-P_{13}^1\left(J_{in}{P}_{in}{T}_{in}{J}_1\left(X_1\right)J_{13}{X}_8{X}_7\right)=0\\ P_{13}^1\left(J_{in}{P}_{in}{T}_{in}{J}_1\left(X_1\right)J_{13}{X}_8{X}_7\right)-P_{21}^1\left(J_{in}{P}_{in}{T}_{in}{J}_2{J}_{21}{X}_9{X}_5\right)=0\\ J_{11}\left(X_2\right)-J_{11}^1\left(X_6\right)=0J_{12}\left(X_3\right)-J_{12}^1\left(X_4\right)=0J_{13}-J_{13}^1\left(X_8\right)=0\\ J_{13}^1\left(X_8\right)-J_{13}^2\left(X_8{X}_7\right)=0J_{21}-J_{21}^1\left(X_5\right)-J_{21}^2\left(X_9\right)=0J_{21}^2\left(X_9\right)-J_{21}^3\left(X_9{X}_5\right)=0\\ 0X_{i}1i=\overline{13}{\epsilon }_{\max }^1{X}_j{\epsilon }_{\max }^{shop}j=\overline{47}{\epsilon }_{\min }^1{X}_k{\epsilon }_{\max }^{1}k=\overline{89}\end{array}$$ E43

where $$P_{11}^1$$ $$P_{12}^1$$ $$P_{13}^1$$ $$P_{21}^1$$ are the natural gas pressure at the outlet of each of the CS branches, $$J_{in}{P}_{in}{T}_{in}$$ are the natural gas flow rate, pressure and temperature at the CS inlet, $$J_1$$ $$J_2$$ are the natural gas mass flow rates through “branches I”, $$J_{11}$$ $$J_{12}$$ $$J_{13}$$ $$J_{21}$$ are the natural gas mass flow rates through “branches II”, $$X_6$$ $$X_4$$ $$X_7$$ $$X_5$$ are the ratios of compression for compressor shops, $$X_8$$ $$X_9$$ are the ratios of compression for GCU groups of the first stages of compressor shops, $${\epsilon }_{\min \max }^1$$ is a minimal/maximal ratio of compression for GCU groups of the first stage, $${\epsilon }_{\max }^{shop}$$ is the maximal ratio of compression for compressor shops.

To assure a safe mode of CS operation, it is required to observe the following restrictions on: maximal volumetric capacity $$Q_j$$ of each operating CFS; frequency $$u_j$$ of the CFS shaft rotation; maximal capacity $$N_j$$ of the CFS drive; maximal outlet pressure $$P_j$$ of the CFS, which is determined by the pipe's strength; maximal temperature $$T_j$$ at the CFS outlet, which is determined by the insulating coating; the lower value of pressure at the outlet of each CFS, related to the requirements to maintain pressure at major gas tapping and boundary gas pipeline points. These restrictions can be formulated as one-sided and two-sided weak inequality:

$$\left\{Q_j^{\min }\le Q_j(\overrightarrow{X})\le Q_j^{\max }{u}_j^{\min }\le u_j\le u_j^{\max }{N}_j(\overrightarrow{X})\le N_j^{\max }{P}_j^{\min }\le P_j(\overrightarrow{X})\le P_j^{\max }{T}_j(\overrightarrow{X})\le T_j^{\max }\right\}$$ E44

where index $j=1\dots M$ means the number of CFS. Also, it is necessary to comply with restrictions on: positions of operational points at the CFS performance curves, related to surge-free operation requirements; conditions related to CFS drive's stable operation; etc. Compliance with the latest restrictions can only be analyzed in the process of a detailed simulation of the mode of gas transmission through the CS.

The system of nonlinear algebraic inequalities (38) describes technological, operating and design constraints on CS equipment operation. When formulating the constraints in a CFD-simulator, a maximum consideration is given to the specifics of the CS process flow chart and the modes of its operation, including a possibility of surging in “CFS – Adjacent CGP” (Seleznev & Pryalov, 2007).

The resultant system (37, 38) represents a generalized system of nonlinear algebraic equalities and inequalities. One of the methods to solve the system of nonlinear algebraic equalities and inequalities is the well-known Interior Point Method, which, in a CFD-simulator, is implemented by statement and solution of an equivalent problem of mathematical optimization (Dennis & Schnabel, 1988). The equivalent optimization problem is solved by the method of modified Lagrange functions. If no solution can be found, the possibility of failure occurrence is admitted. The method we use makes it possible to identify constraints which are not complied with and to evaluate the extent of such non-compliance.

4. Simulation of a GTN by CFD-simulator

Industrial application of CFD-simulators in the gas industry is shown by the example of two acute problems: cutting down the power costs for natural gas transportation through pipeline networks; finding out the reasons for discrepancy between estimates of supplied gas volume by seller and consumer. To resolve these problems, one should be able to simulate the processes of GTN operation in a credible way. Let us consider the GTN simulation algorithm of S. Pryalov as applied to a fragment of a hypothetical CS (fig. 3).

In this CS fragment, three GCUs are combined in parallel into a group by means of a CGP. Fig. 3a uses the following notations: $PT$ – pressure and temperature of the gas (the double number subscripts indicate that the physical parameters of the natural gas flow are attributed to the CFS inlet or outlet with a respective number; parameters at the CS inlet are denoted with the lower index $in$ $J_{i}i=123$ – gas mass flow in the CS pipeline branch with its own number; $n_{i}i=123$ – speed of the $i$ -th CFS shaft; $\xi$ – friction factor for CGP (the double number subscripts indicate that the quantity is attributed to a certain CGP).

Let us conventionally isolate one (the second from the top in fig. 3a) CS branch (fig. 3b): between points $A$ and $B$ . In this case, it is assumed that respective boundary conditions are defined at the branch inlet and outlet points. Such boundary conditions can include: pressure $$P_{\text{A}}$$ and temperature $$T_{\text{A}}$$ at the branch inlet, and pressure $$P_{\text{B}}$$ and heat flux $$q_{\text{B}}$$ (or temperature $T_B$ ) at the branch outlet (note that generally speaking these quantities are functions of time).

Figure 3.

A diagram of a hypothetical CS (a) and a diagram for the branch point A to point B (b)

The value of the mass flow rate through the CFS $$J=J_2^{\text{CFS}}=J_2^{\text{CFS}}\left(t\right)$$ in this case is a solution of a non-linear algebraic equation (lower indices of pressures are the same as in fig. 3):

$${\Psi }_2\left(J\right)=P_{22}^{\text{CFS}}\left(J\right)-P_{22}^{\text{CGP}}\left(J\right)=0$$ E45

where $P_{22}^{\text{CFS}}\left(J\right)={\epsilon }_2\left(P_{21}^{\text{CGS}}\left(J\right)T_{21}^{\text{CGS}}\left(J\right)J\right)\cdot P_{21}^{\text{CGS}}\left(J\right)$ is pressure at the CFS outlet; ${\epsilon }_2\left(P_{21}^{\text{CGS}}\left(J\right)T_{21}^{\text{CGS}}\left(J\right)J\right)$ is gas compression rate produced by the CFS and computed in the CFD-simulator using a semi-empirical CFS model (Seleznev & Pryalov, 2007); $P_{21}^{\text{CGS}}\left(J\right)$ and $T_{21}^{\text{CGS}}\left(J\right)$ are the outlet pressure and temperature of the inlet CGP from the simulations of the gas flow parameters in the given pipeline, in which boundary conditions at its outlet include the target mass flow $J$ and heat flux conservation at the outlet of the inlet CGP; $P_{22}^{\text{CGS}}\left(J\right)=P_{22}^{\text{CGS}}\left(J{T}_{22}^{\text{CFS}}\left(J\right)\right)$ is inlet pressure of the inlet CGP from the simulations of the gas flow parameters in the given pipeline, in which boundary conditions at its outlet include the target mass flow $J$ and gas temperature equal to the gas temperature $T_{22}^{\text{CFS}}\left(J\right)$ at the CFS outlet.

The functions in this equation are determined numerically by using the CFD-simulator at each step of the iterative solution of equation (39). Therefore, the target distributions of the fluid dynamics parameters of both CGPs and operating parameters of the GCU are obtained automatically in the course of solving equation (39).

Equation (39) in the CFD-simulator is solved numerically using the well-known modified Newton procedure, in which equation (39) is solved by iterations using a finite-difference Jacobian approximation. In the paper (Dennis & Schnabel, 1988) the authors show that this method preserves its q-quadratic convergence of the classical Newton procedure for solving non-linear algebraic equations provided that the finite-difference step length is chosen properly.

For simulations of transient CS operation conditions, the computational model of a hypothetical CS fragment (see fig.3a) will comprise sub-models of an inlet CGP, three CFSs and an outlet CGP. The inlet CGP includes pipelines from the point of transported gas flow entry into the CS to the points of the gas flow entry into each CFS. The outlet CGP includes pipelines from the points of gas flow exit from each CFS to the point of gas flow exit from the CS. These inlet and outlet CGPs in this case are branched, multi-line pipelines, which are simulated by models (1–12) or (1, 3, 7, 11, 13–24). It is assumed that boundary conditions at the transported gas flow inlet/outlet points of the CS are defined in accordance with the principles described above. Here, the values of the mass flow rate through each CFS, $$J_i^{\text{CFS}}=J_i^{\text{CFS}}\left(t\right)$$ where $$i=\overline{1N}$$ $$N$$ is the number of CFSs in the group), constitute a solution to the set of non-linear algebraic equations:

$${\Psi }_i\left(J_i\right)=P_{i2}^{\text{CGP}}\left(J_i\right)-P_{i2}^{\text{CFS}}\left(J_i\right)=0i=\overline{1N}$$ E46

The notations used in (40) are similar to those used in equation (39). The subscripts of pressures correspond to the indices shown in fig. 3.

The target distributions of the fluid dynamics parameters of the CGPs and operating parameters of the GCU equipment are obtained automatically by numerical simulation of equations (40) (as a result of using the CFD-simulator for solving this set of equations to determine the values of its constituent functions).

Numerical solution of equations (40) in industrial application of the CFD-simulator is performed by known quasi-Newton methods. In the first line, one can recommend using the modified Broyden method (Dennis & Schnabel, 1988) as one of the best performing extensions of the classical secant method to the case of numerical solution of non-linear algebraic equations. This method can have q-superlinear local convergence with an r-order of $2^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2N$}\right.}$

It is not difficult to extend the above method to the whole GTN. In this case, equations (40) instead of the mass flow rates through the CFSs will contain the mass flow rates through the CS as a whole. By analogy with simulations of gas transportation through an individual CS, due to the use of the CFD-simulator, target distributions of parameters for the whole multi-line GPS and operating parameters of GCU, CFS, DC and ACD are obtained automatically by numerical simulation of modified equations (40) by the quasi-Newton methods.

5. Optimization of gas transmission expenditures using CFD-simulator

To avoid too many technical details without loss of generality, here we consider GTNs with serially connected CSs (fig. 4).

Figure 4.

A diagram of a hypothetical GTN of serially connected compressor stations

Fig. 4 uses the following notations: $PTQ$ – pressure, temperature and volumetric flow rate. The symbols CS1, CS2 and CS3 in fig. 4 are used in the diagram to designate CSs with their serial number. The double number subscripts indicate that the physical parameters of the natural gas flow are associated with the inlet or outlet of the CS with a respective serial number. The horizontal lines between the CSs in the diagram show multi-line GPSs. The flow direction in fig. 4 is from left to right. The GTN diagram is typical of many gas transportation companies. As an example one can mention GTNs of Tomsktransgas, Russia, and SPP, a.s., Slovakia.

The energy criterion of reduced expenses for failure-free natural gas transportation through GTNs has the following form: One should obtain calculated estimates of GTN equipment controls corresponding to minimum costs of integral energy costs related to transient conditions of gas transportation through the GTN in a given time interval subject to compliance with gas specifications at GTN control points and, at the same time, observance of applicable restrictions to ensure safe and environmentally sound operation of the GTN.

In such a criterion, estimates of the cost function for integral energy expenditures are made using the CFD-simulator. For this purpose, the CFD-simulator features functional dependencies of energy carrier consumption in GCUs on the power consumption at CFS shafts for all GCUs of the GTN. Note that plotting distributions of energy carrier consumption in GCUs over time in the CFD-simulator belongs to GTN performance monitoring procedures and is one of its standard functions (Seleznev & Pryalov, 2007). If the energy carrier consumption and the price of these energy carriers are known, it will be easy for a gas transportation company to calculate their total cost.

The requirement of compliance with gas specifications at GTN control points within the criterion is consistent with the requirement for gas transportation companies to fulfill their contractual obligations. The requirement of compliance with restrictions to ensure GTN safety and environmental soundness comes to the observance of process, performance, design and other restrictions as applied to the operation of the GTN equipment. These restrictions can be formalized mathematically as a set of simple constraints on the target values of the GTN equipment controls and a set of non-linear weak two-sided (and/or one-sided) algebraic inequalities. Function values in these inequalities (GTN equipment operating data and fluid dynamics parameters of the gas flow in the network) are calculated using the CFD-simulator.

In practice, in implementing the criterion mathematically, one usually seeks a local minimum of the function of integral expenses subject to given constraints. This is primarily related to the complexity (and, as a rule, polymodality) of the target function, complexity of the constraint functions and uncertainties in input data coming to the CFD-simulator from SCADA systems.

An optimization algorithm for transient GTN operation conditions was developed by V. Seleznev in 2003. It provides for preliminary predictive simulations, using the CFD-simulator, of physical parameters of the gas flow through the GTN for a given time interval based on the predefined history of GTN operation and gas flow parameters at the GTN boundaries. The value of the function of energy (or financial) expenses for gas transportation through the GTN, which is related immediately to the operating parameters of its equipment and gas flow parameters, is calculated at each time step of this simulation. Thus, during the predictive simulations, the time dependence of expenses is automatically generated for a so-called non-optimized prediction.

Setting up an optimized prediction for the given time interval in a general case consists in the definition of time laws for the GTN equipment controls to reduce the costs of natural gas transportation. Target controls should first of all include regulation of the CFS shaft speed in the whole GTN. The speed regulation for all CFSs can result in shutdowns or starts of individual GCUs. The list of involved GCUs, including valve status data, are called the CS configuration (or configuration of the GTN as a whole).

As has been already mentioned above, energy (or financial) expenses at a given time point $$t$$ depend on the status $$\overrightarrow{u}\left(t\right)$$ of the whole package of GTN gas pumping equipment. These expenses will be represented in the form of a non-linear algebraic cost function, $$Z\left(t\overrightarrow{u}\left(t\right)\right)$$ . It is evident that total expenses for the prediction period $$\left[0T\right]$$ will be an integral, $${\displaystyle \underset{0}{\overset{T}{\int }}Z\left(t\overrightarrow{u}\left(t\right)\right)dt}$$ . From the mathematical viewpoint, minimizing the costs of natural gas transportation through the GTN over the predicted period comes to the minimization of the area of the 2D figure $\Phi$ shown in fig. 5.

Figure 5.

Initial function of gas transportation expenses in the analyzed interval (prediction period)

Mathematically, this can be represented in the following way:

$$J\left(\overrightarrow{u}\right)={\displaystyle \underset{0}{\overset{T}{\int }}Z\left(t\overrightarrow{u}\left(t\right)\right)dt}\to \min t\in \left[0T\right]\overrightarrow{u}\left(t\right)\in U\left(t\right)\subset R^m$$ E47

given that

$$\begin{array}{l}\overrightarrow{W}\left(t\overrightarrow{u}\left(t\right)\right)\in \Omega =\{\overrightarrow{W}\in R^k:{\left(w_{\min }^{specified}\right)}_j\le w_j\left(t\overrightarrow{u}\left(t\right)\right)\le {\left(w_{\max }^{specified}\right)}_{j}j=\overline{1l}\\ w_j\left(t\overrightarrow{u}\left(t\right)\right)\le {\left(w_{\max }^{specified}\right)}_{j}j=\overline{l+1k}\}\end{array}$$ E48

where $U\left(t\right)$ is a given set to define the type of constraints related to the limited availability of control resources (simple constraints on the control quantities); $m$ is the number of controls used for optimization; $$\overrightarrow{W}\left(t\overrightarrow{u}\left(t\right)\right)={\left(w_1\left(t\overrightarrow{u}\left(t\right)\right)\mathrm{...}{w}_k\left(t\overrightarrow{u}\left(t\right)\right)\right)}^T$$ is a given algebraic constraint vector-function; $k$ is the total number of constraints in the form of two-sided and one-sided weak algebraic inequalities used for optimization; $${\overrightarrow{w}}_{\min }^{specified}\in R^l$$ and $${\overrightarrow{w}}_{\max }^{specified}\in R^k$$ are given vectors of the minimum and maximum permissible values to implement the constraints; $l$ is the number of constraints in the form of two-sided weak algebraic inequalities used for optimization; $R^m$ is the Euclidean $m$ -dimensional vector space; $R^k$ is the Euclidean $k$ -dimensional vector space; $R^l$ is the Euclidean $l$ -dimensional vector space.

The first $l$ rows of the vector-function $$\overrightarrow{W}\left(t\overrightarrow{u}\left(t\right)\right)$$ describe process, design, environmental and other constraints on the GTN operating conditions. The rest of the constraint vector-function rows provide stabilization of the programmed status and GTN operating conditions. Stabilization of a given programmed status and GTN operating conditions here means maintaining true paths of GTN operating and gas flow parameters in some neighborhood of expected paths. Examples of such GTN parameters are gas pressures at the GTN boundaries, which should be strictly controlled in accordance with contractual obligations of the gas company. Thus, problem (41, 42) includes additional conditions, $$w_j\left(t\overrightarrow{u}\left(t\right)\right)\le {\left(w_{\max }^{specified}\right)}_j$$ $$j=\overline{l+1k}$$ , which will constrain the deviation of computed paths of specified parameters from their programmed paths. Algebraic functions in these constraints are generally modules of discrepancies between computed and programmed paths.

The technology of natural gas transportation through GTNs provides for stepwise variation of the GTN equipment controls (for example, individual GCUs or their groups) in time. Therefore, in setting up an optimal prediction of expenses, one should map a grid $$t_0{t}_1\dots t_N$$ onto the analyzed interval $$\left[0T\right]$$ (see fig. 5) in such a way as to provide constant control of the GTN equipment within the time interval $$\left(t_i{t}_{i+1}\right)$$ $$i=\overline{0N-1}$$ , (fig. 6a), i.e: $$\overrightarrow{u}\left(t\right)={\overrightarrow{u}}^i=const$$ $$t\in \left(t_i{t}_{i+1}\right)$$ $${\overrightarrow{u}}^i\in U\left(t\right)$$ . To simplify the problem definition (without loss of generality), we assume conventionally that the step $$\Delta t=t_{i+1}-t_i$$ $$i=0\mathrm{...}N-1$$ , is constant all over the prediction period.

Figure 6.

Approximation sequence for the cost function

Given this, the function of the GTN equipment status will be represented as follows (fig.6a):

$$u_j(t)=\{\begin{array}{l}\frac{t-t_0}{\delta }\left(u_j^0-u_j\left(0\right)\right)+u_j\left(0\right)t\in \left[0\delta \right)\\ u_j^{i}t\in \left[t_i+\delta t_{i+1}-\delta \right)\\ \frac{t-\left(t_{i+1}-\delta \right)}{2\delta }\left(u_j^{i+1}-u_j^i\right)+u_j^{i}t\in \left[t_{i+1}-\delta t_{i+1}+\delta \right)\end{array}$$ E49

where $$i=0\mathrm{...}N-1$$ $$j=1\mathrm{...}m$$ $$u_j\left(0\right)$$ is a known value of the initial control at the initial time, $$u_j^N=u_j^{N-1}$$ ; the spacing $$\delta$$ is chosen empirically. In industrial applications the function of the GTN equipment status is often a function of shaft speed variation for all CFSs.

It is difficult to calculate the integral of expenses (41) in the constructed optimization model (41–43) analytically. Therefore, this integral is calculated numerically. The widely known method of trapezoids is considered quite acceptable as applied to this case, i.e. (fig.6b):

$$J\left(\overrightarrow{u}\right)={\displaystyle \underset{0}{\overset{T}{\int }}Z\left(t\overrightarrow{u}\left(t\right)\right)dt}\approx 0.5\Delta t{Z}_0\left[\overrightarrow{u}\left(t_0\right)\right]+\Delta t{\displaystyle \sum _{i=1}^{N-1}{Z}_i\left[\overrightarrow{u}\left(t_0+i\cdot \Delta t\right)\right]}+0.5\Delta t{Z}_N\left[\overrightarrow{u}\left(t_N\right)\right]$$ E50

Thus, for some set of discrete controls $${\left({\overrightarrow{u}}^0{\overrightarrow{u}}^1\dots {\overrightarrow{u}}^{N-1}\right)}^T$$ , formula (44) allows calculating the quality functional $$J\left({\overrightarrow{u}}^0{\overrightarrow{u}}^1\dots {\overrightarrow{u}}^{N-1}\right)$$

Hence, the problem of finding an optimal control (41–43) can be replaced with an equivalent problem of nonlinear programming with respect to independent vectors $${\left({\overrightarrow{u}}^0{\overrightarrow{u}}^1\dots {\overrightarrow{u}}^{N-1}\right)}^T$$

$$\tilde{J}\left({\overrightarrow{u}}^0{\overrightarrow{u}}^1\dots {\overrightarrow{u}}^{N-1}\right)=0.5Z_0\left[{\overrightarrow{u}}^0\right]+{\displaystyle \sum _{i=1}^{N-1}{Z}_i\left[{\overrightarrow{u}}^i\right]+0.5Z_N\left[{\overrightarrow{u}}^{N-1}\right]}\to \min$$ E51

subject to the conditions (considering (43)):

$$\begin{array}{l}{\overrightarrow{u}}^i\in U=\left\{\overrightarrow{u}\in R^m:\max {\left[u^{i-1}-{\eta }_{constr}{u}_{\min }^{specified}\right]}_j\le u_j^i\le \min {\left[u^{i-1}+{\eta }_{constr}{u}_{\max }^{specified}\right]}_{j}j=\overline{1m}\right\}\\ i=\overline{1N-1}\end{array}$$ E52

$$\begin{array}{l}{\overrightarrow{W}}_i\left({\overrightarrow{u}}^0{\overrightarrow{u}}^1\dots {\overrightarrow{u}}^{N-1}\right)\in \Omega =\{{\overrightarrow{W}}_i\in R^k:{\left(w_{\min }^{specified}\right)}_j\le w_j^i\left({\overrightarrow{u}}^0{\overrightarrow{u}}^1\dots {\overrightarrow{u}}^{N-1}\right)\le {\left(w_{\max }^{specified}\right)}_j\\ j=\overline{1l}{w}_j^i\left({\overrightarrow{u}}^0{\overrightarrow{u}}^1\dots {\overrightarrow{u}}^{N-1}\right)\le {\left(w_{\max }^{specified}\right)}_{j}j=\overline{l+1k}\}i=\overline{1N}\end{array}$$ E53

where ${\overrightarrow{u}}_{\min }^{specified}\in {\text{R}}^m$ and ${\overrightarrow{u}}_{\max }^{specified}\in {\text{R}}^m$ are given vectors of the minimum and maximum values of control parameters (as a rule, these vectors are time-independent); $${\overrightarrow{\eta }}_{constr}\in {\text{R}}^m$$ is a given constraint on the maximum step of a single j-th control change in time. The numerous constrains in (47) are attributed to the requirement that variation of GTN operating parameters on transition from one time layer to the next one should be smooth. One should be reminded here that calculations of expenses $$Z_i\left({\overrightarrow{u}}^i\right)$$ for gas transportation through the GTN are enabled by the CFD-simulator.

In industrial applications for optimal predictions, knowing the optimal control path is more critical than knowing the minimum value of the quality functional. Therefore, it is reasonable to transform the problem statement (45) in the following way:

$$\stackrel{⌢}{J}\left({\overrightarrow{u}}^1\dots {\overrightarrow{u}}^{N-1}\right)={\displaystyle \sum _{i=1}^N{Z}_i\left({\overrightarrow{u}}^i\right)\to \min }\text{where}{\overrightarrow{u}}^N={\overrightarrow{u}}^{N-1}$$ E54

To solve the problem of general non-linear programming (46–48), we use the method of modified Lagrange functions (Vasilyev, 2002). In accordance with this method, the modified target function of the optimized problem being solved is a sum of initial target function (48) and formalized constrains (47) weighted in a special way (by means of Lagrange multipliers and penalties). An equivalent task of searching for the minimum value of the modified target function subject to simple constraints on variables (46) is fulfilled using modified varied metrics algorithms (Vasilyev, 2002) or modified conjugate gradient algorithms (Vasilyev, 2002), which are resistant to the error accumulation in the course of arithmetic operations.

If the choice of the GTN configuration is considered as an additional control, one will have to introduce additional constraints into the optimization problem to represent process bans on frequent changes of GTN configurations during gas transportation. The choice of the GTN configuration does not require fundamental changes in the optimization approach or algorithm for transient conditions of gas transportation through the GTN.

The above approach to the construction of optimal prediction for gas transportation through the GTN can be easily extended to the CS or group of GCUs.

6. Numerical monitoring of gas distribution discrepancy using CFD-simulator

Numerical monitoring of the discrepancy is based on a statement (for a specified time gap) and numerical solution of identification problem of a physically proved quasi-steady gas dynamics mode of natural gas transmission through specified gas distribution networks.

In large communities, natural gas is supplied to the consumers using medium or low pressure ring mains, being several dozen kilometers long. Gas from the supplier is transmitted to such mains through a GTN after its pressure is reduced by means of a system of gas reducers installed at inlet gas distribution stations (GDSs) (fig. 7). Major parameters of gas supplied by the gas transportation company to the seller are also measured at the GDS outlets. Here, major parameters of natural gas include its flow rate, pressure and temperature.

Gas from inlet GDSs is delivered to the ring main via the CGP network of the gas seller. Consumers receive gas from the ring mains through outlet CGPs leading from the ring main to the consumer. The length of the CGPs can range from several hundred meters to several kilometers. In the first approximation, each consumer is considered independent and provided with gas through one CGP, which is completely associated with the consumer (called “associated CGP” as the text goes). Consumer independence means that the consumer’s gas cannot be delivered to other consumers.

Thus, the gas distribution network (GDN) under consideration comprises inlet CGPs from inlet GDSs, a ring main and associated CGPs. In fig. 7, the GDN under consideration is shown with gray color.

If the GDN operates properly, the seller seeks to sell the whole amount of gas received from the supplier. An exception in this case is natural gas forcedly accumulated in the GDN.

Figure 7.

A diagram of the GDN under consideration

For settlement of accounts, consumers submit reports to the seller, in which they indicate estimated volumes of received gas. These reports are usually generated either by processing the consumers’ field flow meter readings or by simplified calculations based on the rates formally established for the given category of consumers. Verification of data provided by the consumers consists in the comparison of their estimates with data obtained by processing the seller’s flow meter readings in compliance with current guidelines.

The central difficulty in such verification is that the amount of field measurements of supplied gas that can be used as a reliable basis is rather limited in the present-day gas industry. Such a situation results in occasional discrepancies (especially during the heating season) in analyzing the volume of natural gas supplied to the consumers. The total discrepancy over a given time period is determined as a difference between two estimates of the gas volume. The first estimate represents the total gas volume actually received during the time period in question as reported by all consumers, and the second estimate, the total volume of natural gas delivered by the supplier to the seller less the gas volume accumulated in the GDN.

One of the most promising ways to resolve the above problem is to use the CFD-simulator. For this purpose, the following problem setup can be used for numerical monitoring of gas distribution discrepancy using CFD-simulator.

Input data: layout chart of the GDN; sensor locations in the GDN, where gas parameters are measured; given time interval of GDN operation; results of field measurements of gas parameters in the GDN in the given time interval; actual (or nameplate) errors of instruments used to measure gas parameters; data on received gas volumes as reported by each consumer for the given time interval.

Target data: (1) physically based gas flow parameters in the GDN in the given time interval having a minimum discrepancy compared to respective field measurement data at identification points and providing the closest possible agreement between calculated flow rate values at the outlet of each associated CGP and corresponding reported values (further as the text goes, this mode will be called “the identified gas flow”); (2) associated CGPs with underreported gas volumes as against the identified gas flow; (3) calculated estimates of discrepancies between gas volumes delivered in the given time interval through each associated CGP as an arithmetic difference between the calculated gas volume corresponding to the identified gas flow and the reported value; (4) calculated estimates of discrepancies between gas volumes delivered in the given time interval through each inlet GDS as an arithmetic difference between the calculated gas volume corresponding to the identified gas flow and the reported value.

Correct simulation of item 1 in the problem statement makes it possible to obtain credible information on physically consistent space-time distributions of flow rates, pressures and temperatures for the gas flow, which is most reasonable for the given time interval with the given field measurement data. Convergence of calculated and reported gas flow values for individual consumers increases the level of objectivity of numerical analysis, as it seeks to maintain the highest possible trust in the data on received gas volumes reported by the consumers.

It follows from the above problem statement that numerical monitoring of gas distribution discrepancy under items 2–4 in the list of target values in essence consists in performing straightforward arithmetic operations with output data of item 1. Therefore, special attention below will be paid to the algorithm of this calculation. This algorithm was proposed by V. Seleznev in 2008. In the first approximation, we consider the process of gas flow through the GDN to be steady-state.

In order to calculate non-isothermal steady-state gas flow parameters in the GDN under consideration, the following boundary conditions of “Type I” need to be specified: pressure, temperature and composition are defined at the outlet of each inlet GDS; mass flow rate and gas temperature are defined at the outlet of each associated CGP.

Using the CFD-simulator with the given boundary conditions and fixed GDN characteristics, one can unambiguously determine physically based spatial distributions of calculated estimates of steady-state GDN operation parameters (Seleznev et al., 2007). Spatial distributions of parameters here mean their distributions along the pipelines.

A diagram of identification locations is generated on the given layout of sensor locations in the GDN. The preferred location of each identification point should correspond to the key requirement: a considerable change in the fluid dynamics conditions of GDN operation should be accompanied by considerable changes in the gas parameters actually measured at this point. The distribution of identification points over the GDN diagram should be as uniform as possible. An identification point can be located both inside the GDN and at its boundaries. At each identification point, different combinations of major gas flow parameters can be measured. These combinations can be varied for every identification point.

The process of finding the identified gas flow comes to the statement and solution of the problem of conditional optimization:

$${\Vert {\overrightarrow{f}}_{\text{calc}}\left(\overrightarrow{X}\right)-{\overrightarrow{f}}_{\text{meas}}^{\text{const}}\Vert }_L\to \underset{\overrightarrow{X}\in \Omega \subset R^n}{\min }\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}$$ E55

where $${\Vert \dots \Vert }_L$$ is the vector norm, the type of which is determined by the value of the parameter $L$ $\left(L=012\right)$ (see below); $${\overrightarrow{f}}_{\text{calc}}\left(\overrightarrow{X}\right)f_{\text{calc}}:R^n\to R^m$$ is the vector-function of calculated estimates of controlled transported gas variables at the identification points in the $m$ -dimensional Euclidean space $$R^m$$ (these calculated estimates are obtained using the CFD-simulator); $${\overrightarrow{f}}_{\text{meas}}^{\text{const}}\in R^m$$ is a given vector of measured values of controlled transported gas variables at the identification points; $m$ is the number of given identification points in the GDN diagram; $$\overrightarrow{X}\in \Omega \subset R^n$$ is the vector of independent controlled variables in the $n$ -dimensional Euclidean space $$R^n$$ ;

$$\overrightarrow{X}\in \Omega =\left\{\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\overrightarrow{X}\in R^n:\overrightarrow{a}\le \overrightarrow{X}\le \overrightarrow{b}\text{}\text{}{\Vert {\overrightarrow{q}}_{\text{calc}}^{\text{GDS}}\left(\overrightarrow{X}\right)-{\overrightarrow{q}}_{\text{meas\_GDS}}^{\text{const}}\Vert }_0\le t_{\text{flow\_}\text{rate}}^{\text{GDS}}\right\}$$ E56

$$\overrightarrow{a}\in R^n\text{and}\overrightarrow{b}\in R^n$$ are correctly defined vectors setting limits in simple constraints for the range of admissible variation of the vector of independent controlled variables (see below); $n$ is the number of independent controlled variables (see below); $${\Vert \dots \Vert }_0$$ is the cubic vector norm (e.g., $${\Vert \overrightarrow{Y}\Vert }_0=\underset{1\le i\le n}{\max }\left|y_i\right|\overrightarrow{Y}\in R^n$$ $${\overrightarrow{q}}_{\text{calc}}^{\text{GDS}}\left(\overrightarrow{X}\right)q_{\text{calc}}^{\text{GDS}}:R^n\to R^l$$ is the vector function of calculated estimates of mass flow rates through inlet GDSs in the $l$ -dimensional Euclidean space $$R^l$$ (these calculated estimates are obtained using the CFD-simulator); $${\overrightarrow{q}}_{\text{meas\_GDS}}^{\text{const}}\in R^l$$ is a given vector of measured mass flow rates at GDS outlets; $l$ is the number of GDSs; $$t_{\text{flow\_}\text{rate}}^{\text{GDS}}=const$$ is a given upper estimate of actual (nameplate) absolute error of flow meters installed at the inlet GDSs. The constraint in the form of a one-sided weak inequality in (50) formalizes the assumption that the probability of gas underdelivery by the supplier is small.

Components $$x_i$$ of the vector of independent controlled variables here mean some boundary conditions of “Type I” specified in the simulations of steady-state fluid dynamics conditions using the CFD-simulator. As practice shows, problem (49, 50) can be solved successfully, if as components of the vector of independent variables one uses an integrated set of mass flow rates at outlet boundaries of associated CGPs $\left(x_{i}i=\overline{1k}\right)$ and pressures at outlet GDSs $\left(x_{i}i=\overline{k+1n}n=k+l\right)$ , where $k$ is the number of associated CGPs

Components $\left(a_i{b}_{i}i=\overline{1k}\right)$ (see (50)) establish the ranges for controlled variables, the size of which is largely attributed to the degree of the seller’s actual trust in a certain consumer. The following conditions should be necessarily observed:

$$a_i+t_{\text{flow\_}\text{rate}}^{\text{cons}}{\left[q_{cons}^{const}\right]}_i{b}_i-t_{\text{flow\_}\text{rate}}^{\text{cons}}i=\overline{1k}{a}_i+t_{\text{flow\_}\text{rate}}^{\text{cons}}{\left[x_0\right]}_i{b}_i-t_{\text{flow\_}\text{rate}}^{\text{cons}}i=\overline{1k}$$ E57

$${\displaystyle \sum _{i=1}^k{\left[x_0\right]}_i}={\displaystyle \sum _{i=1}^k\left\{{\left[q_{\text{cons}}^{\text{const}}\right]}_i+{\left[\Delta q_{\text{sell}}^{\text{const}}\right]}_i\right\}}={\displaystyle \sum _{j=1}^l{\left[q_{\text{meas\_GDS}}^{\text{const}}\right]}_j}$$ E58

where $${\overrightarrow{q}}_{\text{cons}}^{\text{const}}\in R^k$$ is a given vector of mass flow rates at outlet boundaries of associated CGPs; $$t_{\text{flow\_}\text{rate}}^{\text{cons}}=const$$ is a given upper estimate of the actual (nameplate) absolute error of flow meters installed at outlet boundaries of associated CGPs; $${\overrightarrow{X}}_0\in R^n$$ is the starting point of the conditional optimization problem; $$\Delta {\overrightarrow{q}}_{\text{sell}}^{\text{const}}\in R^k$$ is the increment vector for reported mass flow rates at outlet boundaries of associated CGPs, which is chosen by the gas seller depending on the degree of trust in a certain consumer. Fulfillment of conditions (51) is a guaranty for the consumers that the discrepancy analysis will necessarily account for their reported values of received gas volumes. Constraints (52) serve to implement quasi-steady-state operating conditions of the pipeline network of interest from the very starting point of the conditional optimization problem. The values of remaining components $\left(a_i{b}_{i}i=\overline{k+1n}\right)$ are generally defined in accordance with conditions:

$$a_i+t_{\text{pressure}}^{\text{GDS}}{\left[x_0\right]}_i{b}_i-t_{\text{pressure}}^{\text{GDS}}i=\overline{k+1n}{\left[x_0\right]}_i={\left[p_{\text{meas\_GDS}}^{\text{const}}\right]}_{i-k}i=\overline{k+1n}$$ E59

where $${\overrightarrow{p}}_{\text{meas\_GDS}}^{\text{const}}\in R^l$$ is a given vector of measured pressures in the GDS; $$t_{\text{pressure}}^{\text{GDS}}$$ is a given upper estimate of the actual (nameplate) absolute error of pressure gauges installed in the GDS. As a result of fulfillment of conditions (51) and (53), the starting point of optimization problem (49, 50) will be the inner point with respect to simple constraints on controlled variables, which by far extends the range of methods that can be used for conditional minimization.

Thus, based on (51–53):

$$\begin{array}{l}{a}_i\left\{\min \left({\left[x_0\right]}_i-t_{\text{flow\_}\text{rate}}^{\text{cons}}{\left[q_{\text{cons}}^{\text{const}}\right]}_i-t_{\text{flow\_}\text{rate}}^{\text{cons}}\right)i=\overline{1k}{\left[p_{\text{meas\_GDS}}^{\text{const}}\right]}_{i-k}-t_{\text{pressure}}^{\text{GDS}}i=\overline{k+1n}\right\}\\ b_i\left\{\max \left({\left[x_0\right]}_i+t_{\text{flow\_}\text{rate}}^{\text{cons}}{\left[q_{\text{cons}}^{\text{const}}\right]}_i+t_{\text{flow\_}\text{rate}}^{\text{cons}}\right)i=\overline{1k}{\left[p_{\text{meas\_GDS}}^{\text{const}}\right]}_{i-k}+t_{\text{pressure}}^{\text{GDS}}i=\overline{k+1n}\right\}\end{array}$$ E60

Problem (49–54) can take different forms depending on the type of the vector norm chosen in (49). For example, if we choose the cubic vector norm (L=0), we come to a discrete minimax problem with constraints in the form of one-sided weak inequalities and simple constraints on independent controlled variables:

$$\underset{1\le i\le m}{\max }\left|{\left[f_{\text{calc}}\left(\overrightarrow{X}\right)\right]}_i-{\left[f_{\text{meas}}^{\text{const}}\right]}_i\right|\to \underset{\overrightarrow{X}\in \Omega \subset R^n}{\min }\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}$$ E61

Solution to (55) provides so-called uniform agreement between calculated estimates of gas flow parameters and their measured values. Choosing the octahedron vector norm (L=1) transforms initial problem (49–54) into a general non-linear programming problem represented in the following way:

$$\{\begin{array}{l}{\displaystyle \sum _{i=1}^m\left|{\left[f_{\text{calc}}\left(\overrightarrow{X}\right)\right]}_i-{\left[f_{изм}^{const}\right]}_i\right|}\to \min \\ \overrightarrow{X}\in {\Omega }^{\ast }=\{\overrightarrow{X}\in R^n:\overrightarrow{a}\le \overrightarrow{X}\le \overrightarrow{b}\text{}\text{}\left|{\left[q_{\text{calc}}^{\text{GDS}}\left(\overrightarrow{X}\right)\right]}_j-{\left[q_{\text{meas\_GDS}}^{\text{const}}\right]}_j\right|-t_{\text{flow\_}\text{rate}}^{\text{GDS}}\le 0j=\overline{1l}\}\end{array}$$ E62

Choosing the Euclidean vector norm (L=2) in (49) results in the statement of a new conditional optimization problem, which is almost equivalent to (50-54, 56):

$$\sqrt{{\displaystyle \sum _{i=1}^m{\left({\left[f_{\text{calc}}\left(\overrightarrow{X}\right)\right]}_i-{\left[f_{\text{meas}}^{\text{const}}\right]}_i\right)}^2}}\to \underset{\overrightarrow{X}\in {\Omega }^{\ast }}{\min }$$ E63

Solution of (50-54, 57) provides root-mean-square agreement between calculated estimates of gas flow parameters and their measured values. It should be stressed here that statement (50-55) is stricter than (50-54, 57).

Problems (50-55), (50-54, 56) and (50-54, 57) can be solved numerically using the method of modified Lagrange functions (Vasilyev, 2002), which is quite suitable for this purpose. Note that in practice the time of numerical solution of (50-54, 57) in most cases is much shorter than the time of numerical solution of problems (50–55) or (50-54, 56).

To choose a certain type of the target function in problem (49, 50), a series of numerical experiments were conducted and more than a hundred applied tasks were simulated. The best (in terms of the accuracy/runtime ratio) results in simulating the identification problem (49, 50) were obtained using target function (57).

Based on the above considerations, in order to provide efficiency and improved accuracy of industrial applications, it is reasonable to propose the following algorithm for finding the identified gas flow in the GDN at the initial stage:

Define the starting point $${\overrightarrow{X}}_0\in R^n$$ in accordance with conditions (52) and (53). Define the vectors $$\overrightarrow{a}and\overrightarrow{b}$$ in simple constraints according to (54).

Solve optimization problem (50-54, 57). Results of its numerical solution become input data in searching for the conditional minimum at Step 4.

Analyze correctness of solution results from Step 2. The correctness criterion in this case is the condition of necessary fulfillment of all constraints in problem (50-54, 57). If this criterion is satisfied, proceed to Step 4. If not, extend the variation range of independent variables with subsequent transition to Step 2, i.e. $$\left(\left\{\overrightarrow{a}\overrightarrow{b}\right\}\Rightarrow \left\{{\overrightarrow{a}}^{\ast }{\overrightarrow{b}}^{\ast }\right\}\right)$$ . Usually, the range extension algorithm used here is heuristic and bases on the experience gained in the course of actual simulations.

Find numerical solution to problem (50–55) from the starting point, obtained at Step 2, by the method of modified Lagrange functions. Execution of Step 4 makes it possible to reduce or completely eliminate individual local peaks in discrepancy between calculated estimates and measured values, which may appear at Step 2.

Analyze correctness of the results obtained at Step 4, i.e. check the necessary fulfillment of all constraints in problem (50–55). If the correctness criterion is not fulfilled, solution of Step 3 is assumed to be the target solution.

The vector of controlled variables corresponding to the optimal solution at Step 5, is designated as $${\overrightarrow{X}}_{\text{init}}$$ , with $${\overrightarrow{X}}_{\text{init}}\in \Phi \subset R^n$$ . The found fluid dynamics conditions of GDN operation is taken as the primary fluid dynamics mode. Its calculated parameters have uniform (i.e. strictest) agreement with respective measured values.

At the final stage of identification, the primary fluid dynamics mode is corrected within the available measured information, in order to minimize possible discrepancies between calculated and reported estimates of gas volumes transmitted through each associated CGP in the given time interval. This stage is legal by nature, because given the limited amount of measured data the gas seller has no right to accuse the consumer a priori of deliberate misrepresentation of reported received gas volumes. This stage consists in the solution of the general nonlinear programming problem:

$$\{\begin{array}{l}{\Vert {\overrightarrow{q}}_{\text{calc}}^{\text{cons}}\left(\overrightarrow{X}\right)-{\overrightarrow{q}}_{\text{cons}}^{\text{const}}\Vert }_2\to \min \\ \overrightarrow{X}\in \Theta =\{\overrightarrow{X}\in R^n:a_i^{\ast \ast }\le x_i\le b_i^{\ast \ast }i=\overline{1k}\\ {\left[p_{\text{init\_GDS}}^{\text{const}}\right]}_i-t_{\text{pressure}}^{\text{GDS}}\le x_i\le {\left[p_{\text{init\_GDS}}^{\text{const}}\right]}_i+t_{\text{pressure}}^{\text{GDS}}i=\overline{k+1n}\\ \left|{\left[f_{\text{calc}}^{\text{ident}}\left(\overrightarrow{X}\right)\right]}_s-{\left[f_{\text{init\_ident}}^{\text{const}}\right]}_s\right|-t_{\text{pressure}}^{\text{ident}}\le 0s=\overline{1h}\\ \text{}\text{}\left|{\left[q_{\text{calc}}^{\text{GDS}}\left(\overrightarrow{X}\right)\right]}_j-{\left[q_{\text{meas\_GDS}}^{\text{const}}\right]}_j\right|-t_{\text{flow\_}\text{rate}}^{\text{GDS}}\le 0j=\overline{1l}\}\end{array}$$ E64

where $${\overrightarrow{q}}_{\text{calc}}^{\text{cons}}\left(\overrightarrow{X}\right)\subset \overrightarrow{X}{\left[q_{\text{calc}}\left(x\right)\right]}_i=x_i\text{}i=\overline{1k}$$ is the vector function of calculated mass flow rates for outlet boundaries of associated CGPs in the $k$ -dimensional Euclidean space $$R^k$$ $${\overrightarrow{p}}_{\text{init\_GDS}}^{\text{const}}\in R^l$$ is a given vector of GDS pressure corresponding to the primary fluid dynamics mode at $${\overrightarrow{X}}_{\text{init}}\in R^n$$ $${\overrightarrow{f}}_{\text{calc}}^{\text{ident}}\left(\overrightarrow{X}\right)f_{\text{calc}}^{\text{ident}}:R^n\to R^h$$ is the vector function of calculated estimates of controlled variables at inner identification points in the $h$ -dimensional Euclidean space $$R^h$$ (these calculated estimates are obtained using the CFD-simulator); $${\overrightarrow{f}}_{\text{init\_ident}}^{\text{const}}\in R^h$$ is a given vector of controlled variables at inner identification points corresponding to the primary fluid dynamics mode at $${\overrightarrow{X}}_{\text{init}}\in R^n$$ $h$ is a given number of inner identification points; $$t_{\text{pressure}}^{\text{ident}}=const$$ is a given upper estimate of the actual (nameplate) absolute error of pressure gauges at inner identification points.

The first group of simple constraints on the controlled variables in (58) is partly redundant. It assures that numerical search for solutions in industrial applications is always performed in the domain of practically significant results. The second group of simple constraints and the second group of one-sided weak inequality constraints in problem (58) account for the imperfectness of corresponding existing instruments in favor of consumers. The first group of one-sided weak inequality constraints in (58) formalizes the demand for the closest possible uniform agreement between calculated estimates and reported volumes of gas received by each consumer.

Problem (58) can be solved using modified Lagrange functions (Vasilyev, 2002). As a starting point here we use $${\overrightarrow{X}}_{\text{init}}\in R^n$$ . The target result of the simulation should necessarily be correct, i.e. it should fulfill all simple constraints and inequality constraints of problem (58). Otherwise, the primary fluid dynamics mode is taken as a solution to (58).

The simulation outcome of optimization problem (58) is the final solution to the problem of finding the identification gas flow in the GDN. The target identified gas flow is completely defined by the vector $${\overrightarrow{X}}_{\text{ident}}\in \Theta \subset R^{n+1}$$ , corresponding to the optimal solution of problem (58), and is characterized by the fulfillment of the following conditions: (1) calculated gas flow parameters at each identification point should be as close as possible to corresponding field measurement data; (2) calculated estimates of gas volumes supplied to the GDN in a given time interval should correspond to the supplier-reported values within actual (nameplate) absolute errors of the flow meters installed in the GDS; (3) calculated estimates of gas volumes received by each consumer in the given time interval should be as uniformly close as possible to the values reported by the consumers.

7. Conclusion

The methods and approaches described in this chapter have been developed during the past fifteen years. These methods demonstrated their efficiency as applied to simulations for Gazprom, Russia, and SPP, a.s., Slovakia (Seleznev et al., 2005). These methods apply to pipeline systems for transportation of liquid and gas-liquid products. They are rather efficient in analysis and prevention of accidents (Aleshin & Seleznev, 2004).

The approach presented in this article for high-accuracy numerical analysis of operating parameters of industrial pipeline networks using CFD-simulators is based on adaptation of the full system of equations of fluid dynamics to conditions of transient, non-isothermal processes of the flow of gas mixtures in actual GTNs. The adaptation applies the rule of minimization of the number and depth of accepted simplifications and assumptions. The high accuracy of analysis of industrial pipeline networks operating parameters is understood here as the most reliable description and prediction of actual processes in a GTN, which are achievable due to the present level of development of mathematical modeling and technical monitoring methods and available computer hardware. Development and operation of CFD-simulators in solving industrial problems of improving safety, efficiency and environmental soundness of pipeline network operation can be regarded as one of the promising trends of industrial application of the state-of-the-art computational mechanics methods.

In the past few years, the proposed approach to the numerical monitoring of natural gas delivery through gas distribution networks has proved to be well-performing in industrial applications for discrepancy analysis of natural gas deliveries to large and medium-size consumers. This approach can be implemented even on standard personal computers.

Future development of the methods discussed will be largely focused on their applications in high accuracy analysis and control of multi-phase fluid transportation.