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 [1–4]. 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.
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 can be depicted as a right prism with base area 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. It should be noted that the summarized volume of the joint is equal to, where, is the volume of an infinitely small section of the pipe adjacent to the junction core (see Figure 2(b)). The prism base area can be represented as follows: , where is the factor depending on the prism base geometry only. Now volume can be determined by the following formula:, which means that . The application of this geometrical model enabled us to approximate compliance with mass, momentum, and energy conservation laws at the pipelines junction.
(a)
(b)
Simplifications and assumptions used to construct the heat conductive model of pipeline junction include the following: 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; the simulations take account of “downwind” heat and mass exchange due to heat conduction and diffusion; in the pipeline junction, the gas mixture instantaneously becomes ideally uniform all over the pipeline junction volume (see Figure 2(b)); effects of gas mixture viscosity in the pipeline junction (inside the volume) can be ignored; there are no heat sources in (inside the volume); 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),(ii)for each junction (boundary sections of pipelines adjoining the junction),(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; 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)) are auxiliary function which is necessary for conjugation of axial direction for theth pipeline (r; see Figure 2(b)); are auxiliary weighting function, estimating contribution of theth 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 of the flow cross-section with area is cooling), is the ambient temperature. To denote the relationship of a value to the pipe numbered by, we use a parenthesized superscript on the left side of the value, for example,. 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 is specific (per unit mass) enthalpy of the mixture; 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 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”); associated CGPs with underreported gas volumes as against the identified gas flow; 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: where 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 (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 -dimensional Euclidean space ; are correctly defined vectors setting limits in simple constraints for the range of admissible variation of the vector of independent controlled variables (see below); is the number of independent controlled variables (see below); is the cubic vector norm (e.g., ); is the vector function of calculated estimates of mass flow rates through inlet GDSs in the -dimensional Euclidean space (these calculated estimates are obtained using the CFD simulator); is a given vector of measured mass flow rates at GDS outlets; is the number of GDSs; 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 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 , where is the number of associated CGPs.
Components (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: where is a given vector of mass flow rates at outlet boundaries of associated CGPs; is a given upper estimate of the actual (nameplate) absolute error of flow meters installed at outlet boundaries of associated CGPs; is the starting point of the conditional optimization problem; 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 are generally defined in accordance with conditions where is a given vector of measured pressures in the GDS; is a given upper estimate of the actual (nameplate) absolute error of pressure gauges installed in the GDS. As a result of fulfillment of conditions (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), 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 , 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:
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 transforms initial problem (3.1)–(3.6) into a general nonlinear programming problem represented in the following way: Choosing the Euclidean vector norm 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):
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 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, . 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 , with . 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 -dimensional Euclidean space ; 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 -dimensional Euclidean space (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 ; is a given number of internal identification points; 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 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 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 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] where is a given numerical factor, [14]; is a given numerical factor corresponding to the linear constraint residual decay rate, [15]. We recommend using the following pair [14, 15]: ; .
The initial vector is chosen as close as possible to the optimal vector . For this purpose, we use available a priori information about the solution. The initial value of the parameter 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 . 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.
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.
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.
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: starting point is shifted to a larger (with respect to the basic) flow rate at CGPs; starting point is shifted to a smaller (with respect to the basic) flow rate at CGPs; 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.
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 (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.
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.
In (3.10), the specified upper estimate of the actual (rated) error of pressure sensors at internal identification points was taken equal to, 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 (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.
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 , corresponding to the optimal solution of problem (3.10), and is characterized by the fulfillment of the following conditions: 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.
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.