About this Journal Submit a Manuscript Table of Contents
Journal of Applied Mathematics
Volume 2010 (2010), Article ID 407648, 23 pages
http://dx.doi.org/10.1155/2010/407648
Research Article

Numerical Monitoring of Natural Gas Distribution Discrepancy Using CFD Simulator

Physical & Technical Center, LLC, P.O. Box 236, Sarov, Nizhny Novgorod Region 607190, Russia

Received 31 March 2010; Revised 21 June 2010; Accepted 10 August 2010

Academic Editor: Roberto Natalini

Copyright © 2010 Vadim E. Seleznev. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The paper describes a new method for numerical monitoring of discrepancies in natural gas supply to consumers, who receive gas from gas distribution loops. This method serves to resolve the vital problem of commercial natural gas accounting under the conditions of deficient field measurements of gas supply volumes. Numerical monitoring makes it possible to obtain computational estimates of actual gas deliveries over given time spans and to estimate their difference from corresponding values reported by gas consumers. Such estimation is performed using a computational fluid dynamics simulator of gas flows in the gas distribution system of interest. Numerical monitoring of the discrepancy is based on a statement and numerical solution of identification problem of a physically proved gas dynamics mode of natural gas transmission through specified gas distribution networks. The identified mode parameters should have a minimum discrepancy with field measurements of gas transport at specified reference points of the simulated pipeline network.

1. Introduction

In the last 20 years, operation of complex gas distribution systems has been associated with an acute problem of credible commercial accounting of natural gas supply under the deficiency of respective field measurements [14]. In large communities, natural gas is supplied to the consumers using medium or low pressure ring mains, being several dozen kilometers long. Gasfrom the supplier is transmitted to such mains through a gas transmission network after its pressure is reduced by means of a system of gas reducers installed at inlet gas distribution stations (GDSs) (Figure 1). 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.

407648.fig.001
Figure 1: A diagram of the gas distribution network under consideration.

Gas from inlet GDSs is delivered to the ring main via the connecting gas pipelines (CGPs) 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 Figure 1, 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 represents the total volume of natural gas delivered by the supplier to the seller which is less than the gas volume accumulated in the GDN.

It is necessary to note that in gas industry field measurements [4] were the basic tools for solving a discrepancy monitoring problem for the last decades. The given methods are effective in case of expensive high-precision measuring equipment and providing that flow rate sensors are being installed at all associated CGPs and inlet GDSs of GDN in question. Up to the present day methods of GDN numerical analysis have been of limited application. They were used for GDN geometrical parameters estimation, GDN ultimate capacity analysis, and gas deliveries problem solution [1, 3, 4]. This paper proposes to apply numerical simulation for revealing discrepancy sources of gas delivered through a GDN in the event of lack of measurement information (some of associated CGPs are not equipped with flow rate sensors) and information with a significant level of errors.

The suggested approach to numerical monitoring of the discrepancy is based on a statement and numerical solution of identification problem of a physically proved gas dynamics mode of natural gas transmission through specified gas distribution networks. Identification problem solution is aimed at restoration of space-time distribution of gas flow parameters (pressure, flow rate, and temperature) along the pipelines which constitute GDN in question, on the basis of the incomplete measurement information. The identified mode parameters should have a minimum discrepancy with field measurements of gas transmission at specified reference points of the simulated pipeline network. As such, to numerical evaluation of space-time distribution of gas flow parameters along GDN, one should use the computational fluid dynamics simulator (CFD simulator) [5].

Numerical estimations of gas flow parameters along GDN obtained from the identification problem solution further can be applied for discrepancy monitoring to replace missing measurement information on associated CGPs. Moreover, the comparison of gas volume values calculated and declared by consumers makes it possible to reveal gas volumes which were not paid by a particular consumer. Such approach to gas discrepancy monitoring has never been applied in the gas industry before.

Section 2 describes a mathematical model variant which is used in CFD simulator computational kernel for obtaining calculated estimates of gas transportation parameters.

Section 3 is devoted to statement and numerical solution of identification problem of a physically proved gas dynamics mode of natural gas transmission through specified gas distribution networks.

Other sections of the paper are dedicated to validation and practical application of the approach suggested to numerical monitoring of gas discrepancy.

2. Numerical Kernel of the CFD Simulator

The kernel of the CFD simulator uses a model of gas transmission in branched graded trunk lines [6, 7]. It can be represented as a combination of two types of models: models of gas flow in long pipes adjoining a junction and models of gas flow in junctions.

Branched trunk lines are long, branched, multisection pipelines. For numerical evaluation of parameters of steady and transient, non-isothermal processes of the gas mixture flow in branched trunk lines, 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 [7]. 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. As such, special attention is given to the adequacy of simulation of pipeline junction nodes where the 3D nature of the gas flow is strongly displayed.

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 Pryalov and 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 (Figure 2) proposed by Seleznev et al.[8]. In this model, volume (0)𝑉 can be depicted as a right prism with base area 𝑆base and height 𝐻 (Figure 2(a)). For the prism lateral surface with linear dimensions(𝑛)𝛿, true is the following relation:(𝑛)𝛿=(𝑛)𝑓/𝐻, where (𝑛)𝑓 is the cross-sectional area of the pipe adjacent to the junction core(0)𝑉. It should be noted that the summarized volume of the joint is equal to𝑉=𝑁𝑛=0(𝑛)𝑉, where(𝑛)𝑉, 𝑛=1,,𝑁, is the volume of an infinitely small section of the pipe adjacent to the junction core(0)𝑉 (see Figure 2(b)). The prism base area can be represented as follows: 𝑆base=𝜍base(1)𝛿2, where 𝜍base is the factor depending on the prism base geometry only. Now volume(0)𝑉 can be determined by the following formula:(0)𝑉=𝐻𝑆base=𝐻𝜍base((1)𝑓/𝐻)2=𝜍base(1)𝑓2/𝐻, which means that lim𝐻(0)𝑉=lim𝐻[𝜍base(1)𝑓2/𝐻]=0. The application of this geometrical model enabled us to approximate compliance with mass, momentum, and energy conservation laws at the pipelines junction.

fig2
Figure 2: A Schematic of a pipeline junction: (a) 3D drawing; (b) plain scheme.

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(0)𝑉 (see Figure 2(b)); (4) effects of gas mixture viscosity in the pipeline junction (inside the volume(0)𝑉) can be ignored; (5) there are no heat sources in (0)𝑉 (inside the volume(0)𝑉); (6) pipeline diameters near the pipeline junction are constant.

Then, the heat conductive fluid dynamics model of a transient, non-isothermal, turbulent flow of a viscous, chemically inert, compressible, multicomponent gas mixture through multiline GPS which consist of pipes of round cross-sections and rigid rough heat conductive walls is represented in the following way [7, 8]:(i)for each pipe (bend or nonbranched segment of a ring collector),𝜕(𝜌𝑓)+𝜕𝜕𝑡(𝜕𝑥𝜌𝑤𝑓)=0,(2.1)𝜕𝜕𝑡𝜌𝑌𝑚𝑓+𝜕𝜕𝑥𝜌𝑌𝑚𝜕𝑤𝑓𝜕𝑥𝜌𝑓𝐷𝑚𝜕𝑌𝑚𝜕𝑥=,𝑚=1,𝑁𝑆1,𝑌𝑁𝑆=1𝑁𝑆1𝑚=1𝑌𝑚,(2.2)𝜕(𝜌𝑤𝑓)+𝜕𝜕𝑡𝜌𝑤2𝑓𝜕𝑥=𝑓𝜕𝑝𝜕𝑥+𝑔𝜌𝜕𝑧1𝜋𝜕𝑥4𝜆𝜌𝑤|𝑤|𝑅,(2.3)𝜕𝑤𝜕𝑡𝜌𝑓𝜀+22+𝜕𝑤𝜕𝑥𝜌𝑤𝑓𝜀+22𝜕=𝜕𝑥(𝑝𝑤𝑓)𝜌𝑤𝑓𝑔𝜕𝑧1𝜕𝑥𝑝𝜕𝑓𝜕𝜕𝑡+𝑄𝑓+𝜕𝑥𝑘𝑓𝜕𝑇𝜕𝑥Φ𝑇,𝑇am+𝜕𝜕𝑥𝜌𝑓𝑁𝑆𝑚=1𝜀𝑚𝐷𝑚𝜕𝑌𝑚.𝜕𝑥(2.4)𝜕𝜌+𝜕𝑡𝑁𝑛=1(𝑛)Θ(𝑛)𝜕(𝜌𝑤)𝜕𝑥=0,(2.5)𝜕𝜌𝑌𝑚+𝜕𝑡𝑁𝑛=1(𝑛)𝜕𝜌𝑌𝑚𝑤𝜕𝑥(𝑛)Θ𝑁𝑛=1(𝑛)𝜕𝜕𝑥𝜌𝐷𝑚𝜕𝑌𝑚𝜕𝑥(𝑛)Θ=0,𝑚=1,𝑁𝑆1,(𝑛)𝑌𝑁𝑆=1𝑁𝑆1𝑚=1(𝑛)𝑌𝑚,(2.6)(𝑛)𝜕(𝜌𝑤)+𝜕𝑡(𝑛)𝜕𝜌𝑤2𝜕𝑥=(𝑛)𝜕𝑝𝜕𝑥𝑔𝜌(𝑛)𝜕𝑧1𝜕𝑥0.25(𝑛)𝜆𝜌𝑤|𝑤|𝑅,𝑛=1,𝑁,(2.7)𝜕(𝜌𝜀)+𝜕𝑡𝑁𝑛=1(𝑛)𝜕(𝜌𝜀𝑤)Θ𝜕𝑥=𝑁𝑛=1(𝑛)𝑝𝜕𝑤Θ𝜕𝑥+0.25𝑁𝑛=1(𝑛)𝜆𝜌|𝑤|3Θ𝑅++𝑄𝑁𝑛=1(𝑛)𝜕𝑘𝜕𝑥𝜕𝑇Θ𝜕𝑥𝑁𝑛=1(𝑛)Φ𝑇,𝑇am𝑓Θ+𝑁𝑁𝑛=1𝑠𝑚=1(𝑛)𝜕𝜕𝑥𝜌𝜀𝑚𝐷𝑚𝜕𝑌𝑚Θ,𝜕𝑥(2.8)(𝑛)𝑇=(𝜉)𝑇,𝜀=(𝑛)𝜀=(𝜉)𝜀,(𝑛)𝜀𝑚=(𝜉)𝜀𝑚,𝜌=(𝑛)𝜌=(𝜉)𝜌,(𝑛)𝑝=(𝜉)𝑝,(𝑛)𝑘=(𝜉)𝑘,(𝑛)𝐷𝑚=(𝜉)𝐷𝑚,𝑌𝑚=(𝑛)𝑌𝑚=(𝜉)𝑌𝑚,(𝑛)𝑧1=(𝜉)𝑧1forany𝑛,𝜉1,𝑁,𝑚1,𝑁𝑆,(2.9)𝑁𝑛=1(𝑛)(𝑤𝑓𝑠)=0;𝑁𝑛=1(𝑛)𝜕𝑇𝜕𝑥𝑓𝑠=0;𝑁𝑛=1(𝑛)𝜕𝑌𝑚𝜕𝑥𝑓𝑠=0,(2.10)(𝑛)𝑠=1𝑖𝑓(0)𝐧(𝑛)𝐧<0,1𝑖𝑓(0)𝐧(𝑛)𝐧>0,(𝑛)Θ=(𝑛)𝑉𝑉=(𝑛)𝑓𝐿(𝑛)𝛾𝑁𝑘=1(𝑘)𝑓𝐿(𝑛)𝛾,0<(𝑛)Θ<1,𝑁𝑛=1(𝑛)Θ=1,(2.11)(ii)for each junction (boundary sections of pipelines adjoining the junction),𝑆𝑝=𝑝mix𝑆,𝜀=𝜀mix𝑆;𝑘=𝑘mix,𝜀𝑚=𝜀𝑚𝑆mix,𝐷𝑚=𝐷𝑚𝑆mix,𝑚=1,𝑁𝑆,𝑇1=𝑇2==𝑇𝑁𝑆=𝑇,(2.12)(iii)equations of state (EOS)𝑓

where𝑡 is the density of the gas mixture;𝑥 is the flow cross-sectional area of pipeline; 𝑤 is time (marching variable);𝑌𝑚 is the spatial coordinate over the pipeline's geometrical axis (spatial variable);𝐷𝑚 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);𝑁𝑆 is a relative mass concentration of the m component of the gas mixture;𝑝 is a binary diffusivity of component m in the residual mixture;𝑔 is the number of components of the homogeneous gas mixture;𝑧1 is the pressure in the gas mixture;𝜋 is a gravitational acceleration modulus;𝜆 is the coordinate of the point on the pipeline's axis, measured, relative to an arbitrary horizontal plane, upright;𝑅=𝑓/𝜋 is the Pythagorean number;𝜀 is the friction coefficient in the Darcy-Weisbach formula;𝑄 is the pipe's internal radius;𝑘 is specific (per unit mass) internal energy of the gas mixture;𝑇 is specific (per unit volume) heat generation rate of sources;𝜀𝑚 is thermal conductivity;𝑚 is the temperature of gas mixture;𝑇𝑚 is specific (per unit mass) internal energy of the𝑚 component;𝑁 is the temperature of the(𝑛)𝑠 component;𝑛 is the number of pipes comprising one junction (see (2.5)–(2.11))e(0)𝐧,(𝑛)𝐢 are auxiliary function which is necessary for conjugation of axial direction for the(𝑛)Θth pipeline (r𝑛; see Figure 2(b)); {𝑆mix} are auxiliary weighting function, estimating contribution of theΦ(𝑇,𝑇am)th pipeline for gas flow in junction;𝜒 is a set of parameters of gas mixture. Function𝑓 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(Φ(𝑇,𝑇am)>0 of the flow cross-section with area𝑇am𝑛 is cooling),(𝑛)𝜌 is the ambient temperature. To denote the relationship of a value to the pipe numbered byΦ(𝑇,𝑇am), we use a parenthesized superscript on the left side of the value, for example,Φ(𝑇,𝑇am). In (2.1)–(2.12), we use physical magnitudes averaged across the pipeline's flow cross-section. The set of (2.1)–(2.12) is supplemented by the boundary conditions and conjugation conditions. As to conjugation conditions it is possible to specify boundary conditions simulating a complete rupture of the pipeline and/or its shut-off, operation of valves, and so forth.

As was stated above, the energy equations (2.4) and (2.8) comprise function 𝜀(𝑝,𝑇)=(𝑝,𝑇)𝑝/𝜌; describing the heat exchange between the environment and natural gas in the course of its pipeline transmission. The space-time distribution of function 𝑑=𝑐𝑝𝑑𝑇+(𝜕/𝜕𝑝)𝑑𝑝, 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 [8].

As a rule, Redlich-Kwang equation [4, 8, 9] is used for thermal EOS in (2.12) and calorific EOS—are well-known thermodynamic relations (e.g., 𝑐𝑝 where (1) is specific (per unit mass) enthalpy of the mixture; (2) is heat capacity at constant pressure [8]).

Simulation of steady processes of gas mixture flow through branched trunk lines is a less complicated task compared to (2.1)–(2.12). These models can be easily derived by simplifying the set of (2.1)–(2.12).

To solve (2.1)–(2.12) numerically, the kernel of the CFD simulator generally uses grid methods [8, 10]. Unfortunately, distribution trunk lines contain a large number of CGPs located extremely irregularly at general collectors. In our case, this results in the necessity of considerable spatial grid refinement and, consequently, in a much longer runtime. To resolve this situation, the kernel of the CFD simulator offers a hybrid modification of the known integro-interpolation method [11].

The solution of (2.1)–(2.12) allows to obtain calculated estimates of space-time distribution of gas transportation along GDN. These parameters are used in the net section for identification problem of a physically proved gas dynamics mode of natural gas transmission through specified gas distribution networks.

3. Numerical Detection of Discrepancy Origins

The following problem setup can be used for numerical monitoring of gas distribution discrepancy using CFD simulator.

3.1. Input Data

The input data are 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.

3.2. Target Data

The target data are (3) 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 (control) 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”); (4) associated CGPs with underreported gas volumes as against the identified gas flow; 𝐟mincalc𝐗𝐟constmeas𝐿subjectto𝐗Ω𝑅𝑛,(3.1) 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; 𝐿 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. 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 BC and fixed GDN characteristics, one can unambiguously determine physically based spatial distributions of calculated estimates of steady-state GDN operation parameters. 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:𝐿,(𝐿=0,1,2) where 𝐟calc(𝐗),𝐟calc𝑅𝑛𝑅𝑚, is the vector norm, the type of which is determined by the value of the parameter 𝑚 (see below); 𝑅𝑚 is the vector function of calculated estimates of controlled transported gas variables at the identification points in the 𝐗-dimensional Euclidean space 𝐟constmeas𝑅𝑚 (these calculated estimates are obtained using the CFD simulator with application of (2.1)–(2.12), where 𝑚 is boundary condition at CFD analysis (see below (3.3)–(3.6) and (3.10)); 𝐗Ω𝑅𝑛 is a given vector of measured values of controlled transported gas variables at the identification points; 𝑛 is the number of given identification points in the GDN diagram; 𝑅𝑛 is the vector of independent controlled variables in the 𝐗Ω=𝐗𝑅𝑛𝐪𝐚𝐗𝐛;GDScalc𝐗𝐪constmeas_GDS0𝜏GDSow_rate;(3.2)-dimensional Euclidean space 𝐚𝑅𝑛and𝐛𝑅𝑛;𝑛0 are correctly defined vectors setting limits in simple constraints for the range of admissible variation of the vector of independent controlled variables (see below); 𝐘0=max1𝑖𝑛|𝑦𝑖|,𝐘𝑅𝑛 is the number of independent controlled variables (see below); 𝐪GDScalc(𝐗),𝐪GDScalc𝑅𝑛𝑅𝑙, is the cubic vector norm (e.g., 𝑙); 𝑅𝑙 is the vector function of calculated estimates of mass flow rates through inlet GDSs in the 𝐪constmeas_GDS𝑅𝑙-dimensional Euclidean space 𝑙 (these calculated estimates are obtained using the CFD simulator); 𝜏GDSow_rate=𝑐𝑜𝑛𝑠𝑡 is a given vector of measured mass flow rates at GDS outlets; 𝑥𝑖 is the number of GDSs; (𝑥𝑖,𝑖=1,𝑘) 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 (3.2) formalizes the assumption that the probability of gas underdelivery by the supplier is small.

Components (𝑥𝑖,𝑖=𝑘+1,𝑛,𝑛=𝑘+𝑙) 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 (3.1)-(3.2) 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 𝑘 and pressures at outlet GDSs (𝑎𝑖,𝑏𝑖,𝑖=1,𝑘), where 𝑎𝑖+𝜏consow_rate<𝑞constcons𝑖<𝑏𝑖𝜏consow_rate,𝑖=𝑎1,𝑘,𝑖+𝜏consow_rate<𝑥0𝑖<𝑏𝑖𝜏consow_rate,𝑖=1,𝑘,(3.3)𝑘𝑖=1𝑥0𝑖=𝑘𝑖=1𝑞constcons𝑖+Δ𝑞constsell𝑖=𝑙𝑗=1𝑞constmeas_GDS𝑗,(3.4) is the number of associated CGPs.

Components 𝐪constcons𝑅𝑘 (see (3.2)) 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:𝜏consow_rate=𝑐𝑜𝑛𝑠𝑡 where 𝐗0𝑅𝑛 is a given vector of mass flow rates at outlet boundaries of associated CGPs; Δ𝐪constsell𝑅𝑘 is a given upper estimate of the actual (nameplate) absolute error of flow meters installed at outlet boundaries of associated CGPs; (𝑎𝑖,𝑏𝑖,𝑖=𝑘+1,𝑛) is the starting point of the conditional optimization problem; 𝑎𝑖+𝜏GDSpressure<𝑥0𝑖<𝑏𝑖𝜏GDSpressure,𝑖=𝑥𝑘+1,𝑛,0𝑖=𝑝constmeas_GDS𝑖𝑘,𝑖=𝑘+1,𝑛,(3.5) 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 (3.3) is a guaranty for the consumers that the discrepancy analysis will necessarily account for their reported values of received gas volumes. Constraints (3.4) 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 𝐩constmeas_GDS𝑅𝑙 are generally defined in accordance with conditions𝜏GDSpressure where 𝑎𝑖<𝑥min0𝑖𝜏consow_rate;𝑞constcons𝑖𝜏consow_rate,𝑖=𝑝1,𝑘;constmeas_GDS𝑖𝑘𝜏GDSpressure,𝑖=,𝑏𝑘+1,𝑛𝑖>𝑥max0𝑖+𝜏consow_rate;𝑞constcons𝑖+𝜏consow_rate,𝑖=𝑝1,𝑘;constmeas_GDS𝑖𝑘+𝜏GDSpressure,𝑖=.𝑘+1,𝑛(3.6) is a given vector of measured pressures in the GDS; (𝐿=0) 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 (3.3) and (3.5), the starting point of optimization problem (3.1)-(3.2) 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 (3.3)–(3.5),minmax1𝑖𝑚|||𝑓calc𝐗𝑖𝑓constmeas𝑖|||subjectto𝐗Ω𝑅𝑛.(3.7) Problems (3.1)–(3.6) can take different forms depending on the type of the vector norm chosen in (3.1). For example, if we choose the cubic vector norm (𝐿=1), we come to a discrete minimax problem [12] with constraints in the form of one-sided weak inequalities and simple constraints on independent controlled variables:min𝑚𝑖=1|||𝑓calc𝐗𝑖𝑓constmeas𝑖|||subjectto𝐗Ω=𝐗𝑅𝑛|||𝑞𝐚𝐗𝐛;GDScalc𝐗𝑗𝑞constmeas_GDS𝑗|||𝑡GDSow_rate0,𝑗=.1,𝑙(3.8)

Solution to (3.7) provides so-called uniform agreement between calculated estimates of gas flow parameters and their measured values [13]. Choosing the octahedron vector norm (𝐿=2) transforms initial problem (3.1)–(3.6) into a general nonlinear programming problem represented in the following way:min𝑚𝑖=1𝑓calc𝐗𝑖𝑓constmeas𝑖2subjectto𝐗Ω=𝐗𝑅𝑛|||𝑞𝐚𝐗𝐛;GDScalc𝐗𝑗𝑞constmeas_GDS𝑗|||𝑡GDSow_rate0,𝑗=.1,𝑙(3.9) Choosing the Euclidean vector norm 𝐗0𝑅𝑛 in (3.1) results in the statement of a new conditional optimization problem, which is almost equivalent to (3.2)–(3.6) and (3.8):𝐚and𝐛

Solution of (3.2)–(3.6) and (3.9) provides root-mean-square agreement between calculated estimates of gas flow parameters and their measured values. It should be stressed here that statement (3.2)–(3.7) is stricter than (3.2)–(3.6) and (3.9).

Problems (3.2)–(3.7), (3.2)–(3.6), (3.8), and (3.2)–(3.6), (3.9) can be solved numerically using the method of modified Lagrange functions [14, 15], which is quite suitable for this purpose. Note that in practice the time of numerical solution of (3.2)–(3.6) and (3.9) in most cases is much shorter than the time of numerical solution of problems (3.2)–(3.7) or (3.2)–(3.6) and (3.8).

To choose a certain type of the target function in problem (3.1)-(3.2), 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 (3.1)-(3.2) were obtained using target function (3.9).

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.

Step 1. Define the starting point 𝐚({𝐚,𝐛}{,𝐛}) in accordance with conditions (3.4) and (3.5). Define the vectors 𝐗init in simple constraints according to (3.6).

Step 2. Solve optimization problem (3.2)–(3.6) and (3.9). Results of its numerical solution become input data in searching for the conditional minimum at Step 4.

Step 3. 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 (3.2)–(3.6) and (3.9). If this criterion is satisfied, proceed to Step 4. If not, extend the variation range of independent variables with subsequent transition to Step 2, that is, 𝐗initΦ𝑅𝑛. Usually, the range extension algorithm used here is heuristic and based on the experience gained in the course of actual simulations.

Step 4. Find numerical solution to problem (3.2)–(3.7) 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.

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

Step 6. The vector of controlled variables corresponding to the optimal solution at Step 5 is designated as 𝐪minconscalc𝐗𝐪constcons2subjectto𝐗Θ=𝐗𝑅𝑛𝑎𝑖𝑥𝑖𝑏𝑖,𝑖=𝑝1,𝑘,constinit_GDS𝑖𝜏GDSpressure𝑥𝑖𝑝constinit_GDS𝑖+𝜏GDSpressure,𝑖=|||𝑓𝑘+1,𝑛,identcalc𝐗𝑠𝑓constinit_ident𝑠|||𝜏identpressure0,𝑠=|||𝑞1,,GDScalc𝐗𝑗𝑞constmeas_GDS𝑗|||𝜏GDSow_rate0,𝑗=,1,𝑙(3.10), with 𝐪conscalc(𝐗)𝐗,[𝑞calc(𝑥)]𝑖=𝑥𝑖,𝑖=1,𝑘,. The found fluid dynamics condition 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:𝑘 where 𝑅𝑘 is the vector function of calculated mass flow rates for outlet boundaries of associated CGPs in the 𝐩constinit_GDS𝑅𝑙-dimensional Euclidean space 𝐗init𝑅𝑛; 𝐟identcalc(𝐗),𝐟identcalc𝑅𝑛𝑅, is a given vector of GDS pressure corresponding to the primary fluid dynamics mode at ; 𝑅 is the vector function of calculated estimates of controlled variables at internal identification points in the 𝐟constinit_ident𝑅-dimensional Euclidean space 𝐗init𝑅𝑛 (these calculated estimates are obtained using the CFD simulator); is a given vector of controlled variables at internal identification points corresponding to the primary fluid dynamics mode at 𝜏identpressure=const; 𝐿𝑐(𝐗,𝝁) is a given number of internal identification points; 𝐿𝑐𝐗,𝝁𝑘𝐗+1=𝜔2𝑐𝑘𝑙𝑗=1max0;𝜇𝑗𝑘+𝑐𝑘𝑔𝑗𝐗2𝜇𝑗𝑘2,(3.11) is a given upper estimate of the actual (nameplate) absolute error of pressure gauges at internal identification points.

The first group of simple constraints on the controlled variables in (3.10) 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 (3.10) account for the imperfectness of corresponding existing instruments in favor of consumers. The first group of one-sided weak inequality constraints in (3.10) formalizes the demand for the closest possible uniform agreement between calculated estimates and reported volumes of gas received by each consumer.

Problem (3.10) can be solved using modified Lagrange functions [14]. In accordance with this method, we derive the modified Lagrange function 𝝁𝑘𝑅𝑙:𝑘 where 𝜔(𝐗)=𝑚𝑖=1([𝑓calc(𝐗)]𝑖[𝑓constmeas]𝑖)2 is the vector of Lagrange multipliers at the 𝑐𝑘th iteration of the modified Lagrange functions method; 𝑘 (see (3.9)); 𝐠(𝐗),𝐠𝑅𝑛𝑅𝑙, is the scalar parameter at the 𝐠𝐗=|||𝑞GDScalc𝐗1constmeas_GDS1|||𝑡GDSow_rate|||𝑞GDScalc𝐗𝑙𝑞constmeas_GDS𝑙|||𝑡GDSow_rate𝑅𝑙.(3.12)th iteration of the modified Lagrange functions method; 𝝁𝑘 is the constraint vector function in the form of one-sided slack inequalities in (3.2), that is,𝑘

For the Lagrange multiplier vector 𝑐𝑘 given at the 𝐗𝑘th iteration and the value of the scalar parameter (𝐚𝐗𝐛), the vector 𝜇𝑗𝑘+1=max0;𝜇𝑗𝑘+𝑐𝑘𝑔𝑗𝐗𝐤,𝑗=𝑐1,𝑚+𝑙,𝑘+1=𝛽𝑐𝑘𝐠𝐗𝑖𝑓𝐤0𝐠𝐗>𝛾𝐤𝟏0,𝑐𝑘𝐠𝐗𝑖𝑓𝐤0𝐠𝐗𝛾𝐤𝟏0,(3.13) is defined as a minimum of function (3.11) with simple constraints on the variables 𝛽 (see (3.10)). The minimization problem for function (3.11)-(3.12) with simple constraints on variables can be solved, for example, by the modified conjugate directions method [12], which is stable with respect to the accumulation of arithmetic errors. Then we calculate [14]4𝛽10 where 𝛾 is a given numerical factor, 0<𝛾<1 [14]; 𝛽=10 is a given numerical factor corresponding to the linear constraint residual decay rate, 𝛾=0,25 [15]. We recommend using the following pair [14, 15]: 𝝁0; 𝝁opt.

The initial vector 𝑐0 is chosen as close as possible to the optimal vector 𝐗init𝑅𝑛. For this purpose, we use available a priori information about the solution. The initial value of the parameter (1) should not be too large in order to avoid making the function minimization problem (3.11)–(3.13) artificially ill-conditioned.

As a starting point here we use (2). The target result of the simulation should necessarily be correct, that is, it should fulfill all simple constraints and inequality constraints of problem (3.10). Otherwise, the primary fluid dynamics mode is taken as a solution to (3.10).

4. Method Validation

A series of test simulations was performed to test the performance and efficiency of the method for the numerical detection of discrepancy origins described in Section 3. Let us acquaint ourselves with the statement and results of one test simulation. In this test, we considered a hypothesized GDN. Its computational diagram is shown in Figure 3. Pipeline parameters in the model GDN are shown in Tables 1 and 2.

tab1
Table 1: Pipe parameters in test.
tab2
Table 2: Pipeline parameters in test.
407648.fig.003
Figure 3: A testing diagram of the hypothesized gas distribution network.

As a consistency test of the method developed, we performed a series of simulations of steady-state gas transport in the hypothesized GDN. The boundary conditions for these simulations are given in Table 3. For the purpose of the test simulation, this kind of transport was taken as the basic identified gas dynamic mode.

tab3
Table 3: Boundary Conditions.

The layout of identification points (or reference points) in the hypothesized GDN is shown in Figure 4. Depending on the type of the identification parameter, these points are marked with circles or triangles. Note that the identification points are located not only at the boundaries of the system of interest but also inside it, at Valve 2, Valve 5, and Valve 8. The identification parameters and their values are listed in Table 4.

tab4
Table 4: Gas flow parameters at identification points.
407648.fig.004
Figure 4: A diagram of identification points placement (triangles indicate flow rate monitoring points, circles indicate pressure monitoring points).

Stage 1 of the test simulation is usually performed to test the efficiency of numerical implementation of the initial gas dynamic mode identification algorithm in the CFD simulator kernel. According to the test conditions, the initial state of the given GDN should differ from the basic identified mode. Such an initial state is defined by a starting point, which is characterized by a set of flow rate values at the outlet boundaries of CGPs and gas pressure at GDS outlets. Based on this, one can distinguish three groups of test problems: (3) starting point is shifted to a larger (with respect to the basic) flow rate at CGPs; 3.00 starting point is shifted to a smaller (with respect to the basic) flow rate at CGPs; 3.00 starting point is shifted to a smaller (with respect to the basic) flow rate for some CGPs, and to a larger flow rate for other CGPs. This simulates existing inaccuracies in the consumer-reported values. In addition, the shift in the pressure values at GDS outlets simulates the inaccuracy of sensor-measured values of respective parameters. The objective of the test is that the simulated gas dynamic mode of the system should return to the basic identified mode.

The simulation was performed using the Alfargus/Mosregiongaz computer analytical system (CAS), developed by the Physical & Technical Center, LLC. (Sarov, Russia), which employs the methods described in Sections 2 and 3 (see above). This CAS was based on the Alfargus/Mosregiongaz CFD simulator.There were a total of more than 400 similar runs of Stage 1 of the test simulation. The worst case is represented in Tables 5 and 6.

tab5
Table 5: The worst results after Stage 1 of the test simulation.
tab6
Table 6: Starting point of Stage 1 of the test simulation.

Stage 2 of the test simulation is aimed at revealing the minimum necessary number of internal identification points in the test problem (see Figure 4). For this purpose, we first found the basic mode without internal identification points. Then, internal identification points were introduced one by one by adding gas pressure monitoring sensors at Valve 2, Valve 5, Valve 6, Valve 8, and so forth.

Test results (see Tables 7 and 8) demonstrated that the best closeness to the basic identified mode is obtained for the scheme with three and more internal identification points. The starting point was the same for the first four tests. In the fifth test, the starting point was changed in order to improve reliability of testing. Note that in the tests, the upper estimate of the actual (rated) error of flow meters installed at GDS outlets was taken equal to [kg/s] (see (3.8)). Therefore, the found solutions to tests 1–5 were internal points with respect to the assumed constraints. This evidences the existence of multiple acceptable solutions within the problem statement. However, the identification solution becomes unique and the closest to the basic mode only if there are three and more internal points.

tab7
Table 7: The results after Stage 2 of the test simulation.
tab8
Table 8: Starting point of Stage 2 of the test simulation.

Stage 3 of the test simulation included testing of the third stage of the gas dynamic mode identification problem simulation (see (3.10)), which provides for the correction of the primary gas dynamic mode to minimize possible discrepancies between simulated and reported estimates of the gas volume supplied through each CGP. Thus, the objective of the third stage of testing was to validate the numerical implementation of the maximum consumer credibility criterion.

For this purpose, the result of the first two stages of the identification problem simulation was taken as the primary gas dynamic mode. The resulting mode is shown in Table 9. Here, flow rates reported by some consumers differ from the basic values: flow rates at CGP 4 and CGP 6 (see Table 10). Following this, problem (3.10) is set up and simulated with the starting point of the primary mode.

tab9
Table 9: Starting point of Stage 3 of the test simulation.
tab10
Table 10: Stated mode of natural gas supply to consumers.

In (3.10), the specified upper estimate of the actual (rated) error of pressure sensors at internal identification points 𝑃scal_𝑝=1569064𝑃𝑎 was taken equal to𝑡GDSow_rate=1,5%, where [%] is the given pressure sensor scale value. Also, the upper estimate of the actual (rated) error of flow meters installed at GDS outlets was taken equal to [kg/s] (see (3.10)).

Simulation results presented in Tables 11 and 12 show that the simulated flow rates at CGP 2 and CGP 4 became as close to the reported values as possible subject to the above constraints on monitoring sensors (see the relative deviation for GDS flow meters and absolute deviation for pressure monitoring sensors at CGP 2, CGP 4, Valve 2, Valve 5, and Valve 8). These results indicate that the criterion of (3.10) is satisfied.

tab11
Table 11: The results after Stage 3 of the test simulation.
tab12
Table 12: Final deviations between numerical estimates of gas flow parameters and the results of their measurements.

Stage 1 of the test simulation included more than 100 tests to check the correctness of the CAS kernel operation. Results of these tests were similar to the above results of the first test.

5. On the Aspects of Practical Applications of the Proposed Methods

Note that as the Lagrange particles method (LPM) modification for the energy equation discussed in Section 2 is not related directly to the finite difference grid employed for solving the continuity and motion equations, this grid has almost no effect on the accuracy of the proposed method. Thus, high-accuracy calculated values of gas temperature are obtained without grid refinement, which significantly speeds up the computations.

Due to the absence of direct connection between the LPM method and the finite difference grid, this method is free of the so-called “scheme” viscosity [10]. As a result, the method makes it possible to obtain solutions without scheme smoothing of temperature fronts, which corresponds to real physical processes. This significantly increases the adequacy of simulations compared to the use of difference schemes for the energy equation.

As for the method of numerical detection of discrepancy origins (see Section 3), one should note that it provides the closest possible approximation to consumer-reported gas volume values to within the measurement error of instruments installed at identification points. The simulation outcome of optimization problem (3.10) 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 (2), corresponding to the optimal solution of problem (3.10), and is characterized by the fulfillment of the following conditions: (3) calculated gas flow parameters at each identification point should be as close as possible to corresponding field measurement data; 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; 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.

6. Example of Production Simulations

The above approach to the numerical monitoring of natural gas distribution discrepancy using the CFD simulator is successfully used in the Russian gas industry. Over the period of 2008–2010, this method of numerical monitoring of the supplier share in gas deliveries has demonstrated its efficiency as applied to the production simulations of Mosregiongaz for the analysis of the mechanisms of discrepancy occurrence in the natural gas supplies through the Moscow Ring Gas Distribution Pipeline System (MRGDPS). Figures 5 and 6 show an example of production simulations by the method at Mosregiongaz. The simulation was performed using the CAS system by Mosregiongaz operators jointly with system developers from the Physical & Technical Center, LLC.

407648.fig.005
Figure 5: Numerical simulation of natural gas supply from GDS-17 to Moscow Ring Gas Distribution Pipeline System (an example of instant of time).
407648.fig.006
Figure 6: Numerical simulation of natural gas transmission along MRGDPS (an example of CAS application in Mosregiongaz Control Room).

Different colors in Figure 5 show shares of gas supply from the given GDS. The influence of the GDS on a specific group of consumers is denoted qualitatively by the color of the associated branch, and quantitatively, by the representation of simulation results in the form of tables indicating gas proportions (in percent) consumed by a specific group of consumers from each GDS, which supplies gas to the MRGDPS.

The simulation was performed using CAS by Mosregiongaz operators. Color arrows in the figure indicate the direction and the intensity of gas flows, and spheres of different diameter and color correspond to the actual contribution of each consumer to the discrepancy in MRGDPS gas supply estimates.

7. Conclusion

The developed method enables high-accuracy numerical analysis of discrepancies in natural gas supplies to consumers and scientifically based search for their origins. This method can be fully computerized based on ordinary computers available to gas industry specialists. This makes the method available to gas distribution company specialists, who have no special knowledge in the area of mathematical modeling. With minor modifications, it can be used for the numerical analysis of different networks of trunk and distribution pipelines at power engineering facilities, which transport commercial gases.

References

  1. K. S. Chapman, “Virtual pipeline system testbed to optimize the U.S. natural gas transmission pipeline system,” DOE Report, US Department of Energy, 2003.
  2. T. Kiuchi, “An implicit method for transient gas flows in pipe networks,” International Journal of Heat and Fluid Flow, vol. 15, no. 5, pp. 378–383, 1994.
  3. E. Sekirnjak, “Practical experiences with various optimization techniques for gas transmission and distribution systems,” in Proceedings of the 28th Annual Meeting of the Pipeline Simulation Interest Group (PSIG '96), pp. 1–16, San Francisco, Calif, USA, October 1996.
  4. TAGIF (Technical Association of Gas Industry of France),, Encyclopedia of Gas Industry, 4th edition, 1990.
  5. V. V. Aleshin and V. E. Seleznev, “Computation technology for safety and risk assessment of gas pipeline systems,” in Proceedings of the Asian International Workshop on Advanced Reliability Modeling (AIWARM '04), pp. 443–450, World Scientific, Hiroshima, Japan, August 2004.
  6. S. N. Pryalov and V. E. Seleznev, “Numerical monitoring natural gas supply in gas distribution pipeline pystem,” News of Russian Academy of Science, vol. 1, pp. 152–159, 2010 (Russian).
  7. V. Seleznev, “Numerical simulation of a gas pipeline network using computational fluid dynamics simulators,” Journal of Zhejiang University, vol. 8, no. 5, pp. 755–765, 2007. View at Publisher · View at Google Scholar
  8. V. E. Seleznev, V. V. Aleshin, G. S. Klishin, and R. I. Il’kaev, Numerical Analysis of Gas Pipelines: Theory, Computer Simulation, and Applications, KomKniga, Moscow, Russia, 2005.
  9. A. H. Tarek, Equations of State and PVT Analysis: Applications for Improved Reservoir Modeling, Gulf, Houston, Tex, USA, 2007.
  10. J. D. Anderson Jr., Computational Fluid Dynamics: The Basics with Applications, McGraw-Hill, New York, NY, USA, 1995.
  11. A. A. Samarskii and P. N. Vabishevich, Computational Heat Transfer, John Wiley & Sons, New York, NY, USA, 1995.
  12. J. Nocedal and S. J. Wright, Numerical Optimization, Springer, New York, NY, USA, 2nd edition, 2006.
  13. J. E. Dennis Jr. and R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Ninlinear Equations, Prentice-Hall, Upper Saddle River, NJ, USA, 1988.
  14. D. P. Bertsekas, Nonlinear Programming, Athena Scientific, Belmont, Mass, USA, 2nd edition, 1999.
  15. F. P. Vasilyev, Methods for Optimization, Factorial Press, Moscow, Russia, 2002.