## 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

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

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

for each of the junction nodes

equation of state (EOS) and additional correlations:

where

As was stated above, the energy equations (4) and (8) comprise function

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

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

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

for each of the junction nodes

equation of state:

where the subscript “

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):

where

To record the parametric class of the difference equations (25–30), we used notations of a non-uniform space-time mesh

To explain the notations, it is expedient to consider an individual computational cell containing the node

where*i* node;

Where there was applied the quadratic approximation

the system used the following notations:

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

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:

where

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

Considering the known thermodynamic relationship

where

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:

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.

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

where

To assure a safe mode of CS operation, it is required to observe the following restrictions on: maximal volumetric capacity

where index

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:

Let us conventionally isolate one (the second from the top in fig. 3a) CS branch (fig. 3b): between points

The value of the mass flow rate through the CFS

where

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,

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

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).

Fig. 4 uses the following notations:

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

Mathematically, this can be represented in the following way:

given that

where

The first

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

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

where

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):

Thus, for some set of discrete controls

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

subject to the conditions (considering (43)):

where*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

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:

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.

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:

where

Components

Components

where

where

Thus, based on (51–53):

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:

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:

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):

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

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.

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

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:

where

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

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

## 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.