Numerical Recovery of Gas Flows in Pipeline Systems

Abstract

Optimal control, prevention and investigation of accidents, and detection of discrepancies in estimated gas supply and distribution volumes are relevant problems of trunkline operation. Efficient dealing with these production tasks is based on the numerical recovery of spacetime distribution of nonisothermal transient flow parameters of transmitted gas mixtures based on full-scale measurements in a substantially limited number of localities spaced considerable distances apart along the gas pipelines. The paper describes a practical method of such recovery by defining and solving a special identification problem. Simulations of product flow parameters in extended branched pipelines, involving calculations of the target function and constraint function for the identification problem of interest, are done in the 1D statement. In conclusion, results of practical application of the method in the gas industry are briefly discussed.

1. Problem Statement

Optimal accident-free control of gas trunklines and distribution pipelines, prevention and investigation of accidents in pipeline systems, and detection and localization of discrepancy sources in estimated gas supply and distribution volumes are relevant problems for pipeline transmission, mechanical engineering, and chemical industry [1–6]. Solution of these production problems requires advanced computer simulation methods. These methods use in-depth numerical analysis of commercial and natural gas mixture flow dynamics in high- and medium-pressure linear and circular networks. The cornerstone of such analysis is adequate recovery of spacetime distribution of nonisothermal transient flow parameters of transmitted gas mixtures based on full-scale measurements in a substantially limited number of localities spaced considerable distances apart along the investigated pipeline system.

In other words, solution of effective and accident-free control of pipeline systems requires information about actual pressure, temperature, and gas flow rate distribution along the length of the pipeline and in time. CFD-methods are used for such space-time distributions constructing. At production problems solution one of main reasonableness criteria of obtained numerical estimates of flow parameters is its good correlation in time with full-scale measurements. Meanwhile these measurements are satisfied for limited number of gas flow cross-sections in pipeline. Current cross-sections separate from each other at considerable distance along the length of the pipelines. Modern methods of gas flow modeling along extended branched pipelines can be a base for the obtaining of practically significant correlation between calculated and measured estimates of gas flow parameters by usage of Dirichlet boundary conditions, which are given at chosen boundaries of pipeline network. So, it is reasonable to solve the problem of numerical recovery of gas flows in gas pipeline system by statement and solution of identification problem of fitting calculated and measured estimates of gas flow parameters. The set of changeable in time Dirichlet boundary conditions are assumed as vector of controlled variables of the mentioned problem. The stated above problem is formulated in the form of nonlinear programming problem. As it is known, in practice the global solution of general nonlinear programming problem is rather difficult. That is why it is necessary to develop method of decision of an option with approximate results to the set of full-scale measurements from solution of set of local minimization problems.

As noted above, analysis can be done for linear, branched, or circular pipeline networks. Pipeline segments in such an analysis are composed of single and/or multiple threads made of pipes with rough rigid walls. This paper assumes that the length and location of network pipelines allow for using the one-dimensional setup for the gas network gas flow recovery problem.

Simulated commercial and natural gases are conventionally treated as homogeneous multicomponent heat-conducting viscous gas mixtures of known composition with specified heat transfer, physical and mechanical properties. Equations of state (EOS) for these mixtures are assumed to be known.

Basic modes of their flow through the pipeline networks are assumed to be transient and nonisothermal. At the same time, it is assumed that actual dynamics of real simulated gas flow permits the use of basic assumptions and allowances of the quasi-steady-state flow change method for the gas flow recovery [7, 8].

For the 1D problem statement, one can assume that full-scale measurements of gas flow parameters are taken at fixed points located both at the boundaries of the pipeline system (boundary points) and along the length of the pipelines (internal points). Boundary points are generally used to take measurements of pressure, temperature, and mass flow rate of gases (considering their composition), and internal points are used to measure gas pressure. Ambient temperature is measured at points spaced apart from each other at considerable distances. Results of such measurements may contain random and systematic errors.

In order to avoid too much technical details in the paper, it is assumed in the first approximation that the analyzed gas network contains no injectors, gas compressor units, dust catchers, or gas pressure regulators, and that operation of valves and/or accidents do not alter its configuration over the time interval of interest Δ𝜏. A detailed description of the methods for modeling nominal, transient, and contingency operation of gas trunklines and distribution networks allowing for the operation of valves, compressor and gas distribution plants can be found, for example, in [3, 9–15].

Using the above background information, we should recover distributions of basic gas parameters (i.e., density, pressure, temperature, flow velocity, and composition) along the length of the pipeline and in time for the given interval Δ𝜏.

2. On Modeling of Gas Mixture Flow in the Pipeline Network

Solving the stated problem provides for the use of mathematical models of transient nonisothermal flow of multicomponent gas mixtures developed based on the rule of minimization of the number and depth of necessary adopted assumptions [15]. This is explained by the fact that the use of excessively simplified gas flow models generally leads to the loss of practically essential credibility of simulation results (especially as applied to a detailed analysis of actual gas flow dynamics and search for discrepancies in estimated volumes of natural gas supply to consumers [16]).

An example of a mathematical model satisfying the above requirement is the 1D mathematical model of a transient nonisothermal turbulent flow of a viscous chemically inert compressible multicomponent heat-conducting gas mixture in a branched graded pipeline of circular variable cross-section with absolutely rigid rough heat-conducing walls [17] (Figure 1):

Figure 1 

Schematic representation of a joint of 𝑁 pipelines.

(i) for each pipe:𝜕𝜌𝑓+𝜕𝜕𝑡𝜌𝑤𝑓𝜕𝜕𝑥=0,𝜌𝑌_𝑚𝑓+𝜕𝜕𝑡𝜌𝑌_𝑚𝑤𝑓−𝜕𝜕𝑥𝜕𝑥𝜌𝑓𝐷_𝑚𝜕𝑌_𝑚𝜕𝑥=0,𝑚=1,𝑁_𝑆−1,𝑌_𝑁𝑆=1−^𝑁𝑆−1_𝑚=1𝑌_𝑚,𝜕𝜌𝑤𝑓+𝜕𝜕𝑡𝜌𝑤²𝑓𝜕𝑥=−𝑓𝜕𝑝𝜕𝑥+𝑔𝜌𝜕𝑧₁−𝜋𝜕𝑥4|𝑤|𝜕𝜆𝜌𝑤𝑅,𝑤𝜕𝑡𝜌𝑓𝜀+²2+𝜕𝑤𝜕𝑥𝜌𝑤𝑓𝜀+²2𝜕=−𝑝𝑤𝑓𝜕𝑥−𝜌𝑤𝑓𝑔𝜕𝑧₁𝜕𝑥−𝑝𝜕𝑓𝜕𝜕𝑡+𝑄𝑓+𝜕𝑥𝐾𝑓𝜕𝑇𝜕𝑥−Φ𝑇,𝑇_am+𝜕𝜕𝑥𝜌𝑓^𝑁𝑆_𝑚=1𝜀_𝑚𝐷_𝑚𝜕𝑌_𝑚;𝜕𝑥 (2.1)(2.2)(2.3)(2.4)

(ii) junction conditions for each joint: ^𝑁(𝑛)_𝑛=1𝑤𝑓_𝐿𝑠=0;^𝑁(𝑛)_𝑛=1𝑓𝜕𝑇𝜕𝑥_𝐿𝑠=0,_{𝑁(𝑛)𝑛=1}𝑓𝜕𝑌_𝑚𝜕𝑥_𝐿𝑠=0,𝑚=1,𝑁_𝑆−1,_(𝑛)𝑇_𝐿=^(𝜉)𝑇_𝐿,𝜀_Joint=^(𝑛)𝜀_𝐿,^(𝑛)𝜀_𝑚𝐿=^(𝜉)𝜀_𝑚𝐿,𝜌_Joint=^(𝑛)𝜌_𝐿,𝑝_Joint=^(𝑛)𝑝_𝐿,_(𝑛)𝑘_𝐿=^(𝜉)𝑘_𝐿,^(𝑛)𝐷_𝑚𝐿=^(𝜉)𝐷_𝑚𝐿,𝑌_𝑚,Joint=^(𝑛)𝑌_𝑚𝐿,^(𝑛)𝑧_1𝐿=^(𝜉)𝑧_1𝐿∀𝑛,𝜉∈1,𝑁and∀𝑚∈1,𝑁_𝑆; (2.5)(2.6)

(iii) fitting conditions for each joint: 𝜕𝜌_Joint+𝜕𝑡^𝑁(𝑛)_𝑛=1Θ^(𝑛)𝜕𝜌𝑤𝜕𝑥_𝐿𝜕=0,𝜌𝑌_𝑚Joint+𝜕𝑡^𝑁(𝑛)_𝑛=1Θ^(𝑛)𝜕𝜌𝑌_𝑚𝑤−𝜕𝜕𝑥𝜕𝑥𝜌𝐷_𝑚𝜕𝑌_𝑚𝜕𝑥_𝐿𝑌=0,_{𝑁𝑆,Joint}=1−^𝑁𝑆−1_𝑚=1𝑌_𝑚,Joint,𝑚=1,𝑁_𝑆−1,_(𝑛)𝜕𝜌𝑤+𝜕𝜕𝑡𝜌𝑤²+𝜕𝑥𝜕𝑝𝜕𝑥+𝑔𝜌𝜕𝑧₁|𝑤|𝑅𝜕𝑥+0.25𝜆𝜌𝑤⁻¹_𝐿=0,𝑛=𝜕1,𝑁,𝜌𝜀_Joint+𝜕𝑡^𝑁(𝑛)_𝑛=1Θ^(𝑛)𝜕𝜌𝜀𝑤𝜕𝑥+𝑝𝜕𝑤|𝑤|𝜕𝑥−0.25𝜆𝜌³𝑅⁻¹−𝜕𝐾𝜕𝑥𝜕𝑇𝜕𝑥_𝐿+^𝑁(𝑛)_𝑛=1Θ^(𝑛)Φ𝑇,𝑇_am𝑓−^𝑁𝑆_𝑚=1𝜕𝜕𝑥𝜌𝜀_𝑚𝐷_𝑚𝜕𝑌_𝑚𝜕𝑥−𝑄_𝐿=0; (2.7)(2.8)(2.9)(2.10)

(iv) EOS and additional relations: 𝑝=𝑝𝑆_mix,𝜀=𝜀𝑆_mix,𝐾=𝐾𝑆_mix,𝑇₁=𝑇₂=⋯=𝑇_𝑁𝑆𝜀=𝑇,_𝑚=𝜀_𝑚𝑆_mix,𝐷_𝑚=𝐷_𝑚𝑆_mix,𝑚=1,𝑁_𝑆; (2.11)

(v) auxiliary geometric relations: ^(𝑛)𝑠=−^(0)→𝐧⋅^(𝑛)→𝐢=1,if^(0)→𝐧⋅^(𝑛)→𝐢<0,−1,if^(0)→𝐧⋅^(𝑛)→𝐢>0,_(𝑛)Θ=^(𝑛)𝑉𝑉=^(𝑛)𝛾𝑓∑^𝑁(𝑘)_𝑘=1𝛾𝑓,0<^(𝑛)Θ<1,_{𝑁(𝑙)𝑙=1}Θ=1;^(𝑛)𝛾=^(𝑛)Δ𝑋Δ𝑋,𝑛=1,𝑁;_(𝑛)Δ𝑋=^(𝑛)𝑥_𝐿−^(𝑛)𝑥_𝑀=^(𝑛)𝑠^(𝑛)𝑥_𝐿−^(𝑛)𝑥_𝑀=^(𝑛)𝑠𝑥_𝐿−𝑥_𝑀. (2.12)(2.13) The function Φ(𝑇,𝑇_am) characterizes the heat exchange of the gas flow core through the boundary gas layer, pipe wall, and insulation with the environment. It expresses a specific (per unit length) total thermal flux along the perimeter 𝜒 of the cross-section having an area of 𝑓 from the transported gas to the environment (Φ(𝑇,𝑇_am)>0 corresponds to heat removal; 𝑇_am is the spacetime distribution of ambient temperature at the domain boundary). The system of equations (2.1)–(2.13) is supplemented with boundary conditions.

The first version of the model (2.1)–(2.13) was proposed by Seleznev et al. jointly with the author of this paper at the end of the last century in order to improve the credibility of gas trunkline flows modeling by suppressing the adverse effects of numerical discrepancies that are completely difference-based [18]. The model (2.1)–(2.13) uses the fitting conditions (2.7)–(2.10), which—together with the junction conditions (2.5) and (2.6)—serve to provide guarantied fulfillment of (2.1)–(2.4) in pipeline joints in their numerical analysis involving grid methods (including simulations on coarse spatial grids). This provides correct—from the practical point of view—compliance with the basic conservation laws in branched pipeline networks, including the mass balance of the transported product. Many years of active use of this model in production simulations of pipeline system confirmed its high efficiency [19–24]. In addition, it should be noted that the model of steady-state nonisothermal gas mixtures transmitted in pipeline network can be obtained through straightforward operations to simplify the system of equations (2.1)–(2.13) using natural assumptions that partial time derivatives are equal to zero.

The model (2.1)–(2.13) implicitly includes the parameters ^(𝑛)𝛾,𝑛=1,𝑁 (2.13). It is evident that these parameters are strictly geometric characteristics. In addition, the choice of their specific values is not constrained in any way. As a result, it may seem that the solution of the system of partial differential equations (2.1)–(2.13) depends on arbitrarily chosen geometric parameters ^(𝑛)𝛾,𝑛=1,𝑁, which makes it ambiguous. However, one can demonstrate that this model remains correct for arbitrary values of^(𝑛)𝛾,𝑛=1,𝑁. Validation of this statement is beyond the scope of this paper, but its detailed description can be found in [17, 25].

A detailed presentation of the practical methods for numerical solution of the system (2.1)–(2.13) and its modification for steady-state flows can be found, for example, in [15, 17, 25–27].

3. Formalization of the Gas Mixture Flow Recovery Problem and Method of Its Solution

The procedure of transmitted flows recovery for linear and/or circular pipeline systems based on “noisy” field measurement data can be formalized as a statement and solution of a special mathematical identification problem. For this purpose, we introduce the notion of identification point (IP). In our case, the IP is an inner or boundary point in the computational model of the pipeline network of interest, in which full-scale measurements of pressure of the transported gas mixture are taken over a given time interval Δ𝜏. Note that the computational model of a gas pipeline system is built according to the principle of minimization of assumptions made in the description of the system’s real topology supported by actual or rated characteristics of its segments. The choice of gas mixture pressure as an identification parameter is explained by the fact that pressure histories in real pipeline systems are determined today more accurately than temperature or flow rate parameters.

In the course of mathematical identification, calculated and measured estimates of gas mixture pressure histories for the entire set of IPs distributed across the computational model of the pipeline network should be fitted as closely as possible. The preferable location of each IP is determined subject to the following requirement: any considerable change in gas dynamic parameters of pipeline system operation should be accompanied by considerable changes in gas mixture flow parameters actually measured at this point. The distribution of IPs in the computational model of the pipeline system should be as uniform as possible. The close fit between corresponding calculated and measured pressure histories in the general case should be provided in three senses [28, 29]: (1) close fit between two functional relations (in essence, between the first derivatives of the functions being compared); (2) close fit between two functional relations in the time-weighted average metric 𝐿₁ (in our case, it is defined by means of the octahedral vector norm, that is, 𝐿₁=‖𝐘^→‖₁=∑^𝑛_𝑖=1|𝑦_𝑖|,𝐘^→∈𝑅^𝑛) or metric 𝐿₂ (in our case, it is defined by means of the Euclidean vector norm,that is, 𝐿₂=‖𝐘^→‖₂=[∑^𝑛_𝑖=1𝑦²_𝑖]^0.5,𝐘^→∈𝑅^𝑛); (3) close fit between two functional relations within their uniform deviation, that is, in the metric 𝐿₀ (in our case, it is defined by means of the cubic vector norm, that is, 𝐿₀=‖𝐘^→‖₀=max_{1≤𝑖≤𝑛}|𝑦_𝑖|,𝐘^→∈𝑅^𝑛).

Real pipeline systems contain a number of branches, through which transmitted fluids enter or leave the system. Inlet branches that supply the gas mixture into the simulated pipeline system will be designated conventionally as “supplier branches,” and outlet branches, as “consumer branches.” In the first approximation, it is assumed that each network branch cannot change its purpose over a given time interval Δ𝜏, that is, a gas supplier cannot become a consumer and vice versa.

In practice, there is generally a shortage of instruments at outlet boundaries of the gas pipeline system of interest. In this case, a number of consumers having no flow rate meters joint declare their gas consumption based on regulatory documents. All the foregoing (together with real instrument errors and encountered cases of artificial under-/overdeclaration of gas mixture volumes transported through the pipeline system) results in arithmetic discrepancies between estimated volumes of gas supply made by consumers and suppliers. This situation should be taken into account in gas flows recovery.

Given all the reasoning above, the special problem of mathematical identification can be stated using conditional optimization: 𝑃𝑡,𝐹_scenario,^→𝐙(𝑡)→min_{→𝐙(𝑡)(𝑡)∈𝚷⊂𝑅𝑛}. (3.1)

In our case, components of the vector-function 𝐙^→(𝑡) are time functions of pressure and mass flow rates of the gas mixture at outlet boundaries of the gas pipeline system, that is, components of boundary conditions. They are distributed in the following way: gas mixture pressures at 𝑢(𝑢≤𝑙) outlets of supplier branches and at 𝑠(𝑠≤𝑘) outlets of consumer branches, where 𝑘 is a given number of consumer branches, through which the gas mixture leaves the gas pipeline system over the time Δ𝜏; mass flow rates of the gas mixture at (𝑙−𝑢) outlets of supplier branches and at (𝑘−𝑠) outlets of consumer branches. Hence, the total number of variable components in the boundary conditions is 𝑛=𝑙+𝑘. In production simulations, the search for boundary conditions at outlet boundaries of a number of branches during identification can be replaced with rigidly set Dirichlet boundary conditions in the form of combinations of known time functions of measured mass flow rates and pressures of the gas mixture. This reduces the number of variable components in the boundary conditions, 𝑛<𝑙+𝑘. Note also that as Dirichlet boundary conditions for temperature and relative mass fraction one can use predefined time laws for respective measured values.

When running production simulations, pipeline system operators face the necessity of gas flow recovery under three basic assumptions: estimated volumes of gas mixture supply declared by suppliers and consumers contain errors (the Full Distrust computational scenario); only supplier-declared estimated volumes of gas mixture supply are credible (the trust-in-supplier computational scenario); only consumer-declared estimated volumes of gas mixture supply are credible (the trust-in-consumer computational scenario). The flag of the computational scenario assigned to the identification problem 𝐹_scenario takes the values of 11, 12, 13, 21, 22, 23, 31, 32, and 33 in series, allowing us to choose various modifications of the problem statement (3.1). The way of using the set of values assigned to the flag 𝐹_scenario will be demonstrated below (see (3.2)–(3.9)).

The target function 𝑃[𝑡,𝐹_scenario,^→𝐙(𝑡)] of the problem (3.1) subject to the requirement that calculated and measured values should fit together in the three above senses can be formalized as follows:𝑃𝑡,𝐹_scenario,^→𝐙(𝑡)=^→𝐩^IP_calc𝑡,𝐙^→(𝑡)−^→𝐩^IP_meas(𝑡)₂,if𝐹_scenario𝜕<21,𝜕𝑡^→𝐩^IP_calc𝑡,^→𝐙(𝑡−𝜕)𝜕𝑡^→𝐩^IP_meas(𝑡)₂,if𝐹_scenario<31,_→𝝎_I-II𝑡,𝐙^→(𝑡)₂,otherwise,_→𝝎_I-II𝑡,𝐙^→(𝑡)=^→𝐩^IP_calc𝑡,^→𝐙(𝑡)−^→𝐩^IP_meas(𝑡)_I+^→𝐩^IP_meas(𝑡)−^→𝐩^IP_calc𝑡,^→𝐙(𝑡)_II. (3.2) The vector-function ^→𝐩^IP_calc[𝑡,^→𝐙(𝑡)] is built by numerical solution of (2.1)–(2.13) for the known initial conditions and defined Dirichlet boundary conditions, containing all the components of the vector-function 𝐙^→(𝑡).

The first form of the target function (i.e., for 𝐹_scenario<21 in (3.2)) in the problem statement (3.1) expresses the requirement that calculated and measured estimates of gas mixture pressure should be close in the second and third senses (see above). In practice, striving for the fulfillment of this requirement makes it possible to obtain a correct solution in the presence of random errors in flow pressure measurements aggravated by single instrument failures. Note that the results of minimization of (3.1), when searching for the sources of discrepancies between natural gas supply estimates with such a form of the target function, are most reasonable from the legal point of view.

The second form of the target function (for 𝐹_scenario<31 in (3.2)) was proposed to enable a closer fit between calculated and measured estimates of gas mixture pressure in the first and third senses (see above). This requirement is aimed at obtaining a correct solution in the presence of systematic errors in gas mixture pressure measurements and single instrument failures. It makes sense to note that under production simulations it failed to fit calculated and measured estimates of gas mixture pressure in all the three senses at once. This was primarily attributed to a fairly high level of “noise” (random and systematic errors) field measurement data.

The third form of the target function in (3.2) was proposed by V. V. Kiselev, first of all, in order to compensate for the shortage of IPs, which is frequently experienced in practice. As noted in the legend to (3.2), both natural (physically based) pressure differences between neighbor IPs (including those in circular pipelines) and virtual pressure differences between intentionally chosen pairs of IPs are considered here. IPs in each pair of points determining the controlled virtual pressure differences can be located in far parts of the analyzed pipeline system. In order to define natural and virtual gas mixture pressure differences controlled during minimization of (3.1), a generalized set of IP pairs is established in advance. The composition of this set does not depend on time. It is established using simple rules: IP pairs should include (at least once) each IP in the computational model of the pipeline system; repetitions and one of “mirror-reflected” IP pairs (if present in the original set) is necessarily excluded from the set; the condition 𝑀_difference≫𝑀_IP should be strictly fulfilled. If the third rule is violated, it seems unreasonable to use the third form of the target function for production simulations.

Now, let us proceed to discussing numerous constraints on (3.1):𝐙^→(𝑡)∈𝚷(𝑡)=𝐙^→(𝑡)∈𝑅^𝑛∶^→𝐠(𝑡)≤^→𝐙(𝑡)≤^→𝐟(𝑡),𝑔_𝑞(𝑡)_𝑠≤𝑞^supplier_calc𝑡,𝐙^→(𝑡)_𝑠≤𝑓_𝑞(𝑡)_𝑠,𝑠=𝑔1,𝑙,_𝑞(𝑡)_𝑠≤𝑞^consumer_calc𝑡,^→𝐙(𝑡)_𝑠−𝑙≤𝑓_𝑞(𝑡)_𝑠,𝑠=𝑙+1,𝑙+𝑘,unequalityΔ𝜏,𝐹_scenario,𝐙^→(𝑡). (3.3) The vector-function ^→𝐪^consumer_calc[𝑡,𝐙^→(𝑡)] is built by numerical solution of (2.1)–(2.13) for the known initial conditions and defined Dirichlet boundary conditions, containing all the components of 𝐙^→(𝑡). The following law of signs is true for the vector-function ^→𝐪^consumer_calc[𝑡,𝐙^→(𝑡)]: if [𝑞^consumer_calc[𝑡,^→𝐙(𝑡)]]_𝑖<0, the gas mixture moves from the gas pipeline system to the consumer, 𝑖=1,𝑘. The second and third inequalities in the list of constraints make it possible to reliably control the variations in gas mass flow rates at outlets of all system branches irrespective of whether these functions are components of the vector-function of controlled boundary conditions, or they are purely computational parameters needed for simulations of gas mixture flow through the pipeline system.

The formal representation of the inequality in (3.3) can be expanded in the following way: unequalityΔ𝜏,𝐹_scenario,^→𝐙(𝑡=)−∫^𝑘_Δ𝜏∑_𝑗=1𝑞^consumer_meas(𝑡)_𝑗𝑑𝑡−Δ≤∫^𝑙_Δ𝜏∑_𝑖=1𝑞^supplier_calc𝑡,^→𝐙(𝑡)_𝑖𝑑𝑡≤∫^𝑙_Δ𝜏∑_𝑖=1𝑞^supplier_meas(𝑡)_𝑖𝑑𝑡+Δ,if∫^𝑙_Δ𝜏∑_𝑖=1𝑞^supplier_meas(𝑡)_𝑖𝑑𝑡≥−∫^𝑘_Δ𝜏∑_𝑗=1𝑞^consumer_meas(𝑡)_𝑗𝑑𝑡,𝐹_scenario,∫=11or21or31_{𝑙Δ𝜏𝑖}∑₌₁𝑞^supplier_meas(𝑡)_𝑖𝑑𝑡−Δ≤∫^𝑙_Δ𝜏∑_𝑖=1𝑞^supplier_calc𝑡,𝐙^→(𝑡)_𝑖𝑑𝑡≤−∫^𝑘_Δ𝜏∑_𝑗=1𝑞^consumer_meas(𝑡)_𝑗𝑑𝑡+Δ,if∫^𝑙_Δ𝜏∑_𝑖=1𝑞^supplier_meas(𝑡)_𝑖𝑑𝑡<−∫^𝑘_Δ𝜏∑_𝑗=1𝑞^consumer_meas(𝑡)_𝑗𝑑𝑡,𝐹_scenario=11or21or31_→,𝐪^consumer_calc𝑡,𝐙^→(𝑡)−^→𝐪^consumer_meas(𝑡)₂−𝜀≤0,if𝐹_scenario=13or23or33,_→𝐪^supplier_calc𝑡,𝐙^→(𝑡)−^→𝐪^supplier_meas(𝑡)₂−𝜀≤0,otherwise. (3.4) The vector-function ^→𝐪^supplier_calc[𝑡,^→𝐙(𝑡)] is built by numerical solution of (2.1)–(2.13) for the known initial conditions and Dirichlet boundary conditions, containing all the components of the vector-function of the controlled boundary conditions. In this case, the following law of signs is valid for the components of the mass flow rate vector-functions introduced above. [𝑞^supplier_calc[𝑡,𝐙^→(𝑡)]]_𝑖>0 or [𝑞^supplier_meas(𝑡)]_𝑖>0 means that the gas mixture enters the gas pipeline system, 𝑖=1,𝑙. The following law of signs is true for components of the vector-function ^→𝐪^consumer_meas(𝑡): if [𝑞^consumer_meas(𝑡)]_𝑖<0, the gas mixture moves from the gas pipeline system to the consumer, 𝑖=1,𝑘.

Today, it does not seem possible to solve the problem (3.1)–(3.4) in such a statement using computing facilities available to a wide range of pipeline industry specialists. However, as mentioned in Section 2, actual operation dynamics of most commercial gas trunklines renders it possible to use basic allowances and assumptions of the quasi-steady-state flow change method. In this connection, it is suggested that the time interval of interest Δ𝜏 be conventionally divided into (𝑁_𝑡+1) time layers separated from each other by a given uniform step Δ𝜏. The 𝑚=0 layer will correspond to the lower boundary of the time interval Δ𝜏, and the 𝑚=𝑁_𝑡 layer, to its upper boundary. In order to improve the credibility of estimated gas mixture supply to consumers, when using the quasi-steady-state (for one time layer) problem statement, one should give consideration to the effect of product buildup in the pipes of the simulated pipeline system. For each time layer, the gas mixture buildup varies over the preceding time interval Δ𝑡. A practical way of accounting for this buildup was proposed by Seleznev et al. [17, 30]. As noted above, quasi-steady-state pipeline system operation parameters can be calculated with a simplified version of the model (2.1)–(2.13).

Thus, using the above assumptions and applying the difference approximation of partial time derivatives present in (3.1)–(3.4), the special identification problem can be stated in the following discrete form: 𝑃𝐹_scenario,^→𝐙_𝑚→min_{→𝐙𝑚∈𝚷𝑚⊂𝑅𝑛};𝑃𝐹_scenario,^→𝐙_𝑚=^𝑀IP_𝑖∑₌₀𝑝^IP_calc^→𝐙−𝑝^IP_meas^𝑖_𝑚−Ω_𝑚𝐹_scenario𝑝^IP_calc^→𝐙−𝑝^IP_meas^𝑖_𝑚−1²,if𝐹_scenario<31;_{𝑀difference𝑖}∑₌₀𝑝^IP_calc^→𝐙−𝑝^IP_meas^I_𝑚^,𝑖+𝑝^IP_meas−𝑝^IP_calc^→𝐙^II_𝑚^,𝑖2,otherwise,_→𝐙_𝑚∈𝚷_𝑚=^→𝐙_𝑚∈𝑅^𝑛∶^→𝐠_𝑚≤𝐙^→_𝑚≤^→𝐟_𝑚,𝑔_𝑞^𝑠_𝑚≤𝑞^supplier_calc^→𝐙^𝑠_𝑚≤𝑓_𝑞^𝑠_𝑚,𝑠=𝑔1,𝑙,_𝑞𝑠𝑚≤𝑞^consumer_calc^→𝐙_𝑚^𝑠−𝑙≤𝑓_𝑞^𝑠_𝑚,𝑠=𝐹𝑙+1,𝑙+𝑘,unequality_scenario,𝐙^→_𝑚,𝑚=0,𝑁_𝑡, (3.5)(3.6)(3.7) where Ω_𝑚𝐹_scenario=0,if𝑚=0or𝐹_scenario𝐹<21,1,otherwise,unequality_scenario,^→𝐙_𝑚=−∑^𝑘_𝑗=1𝑞^consumer_meas^𝑗_𝑚−Δ≤∑^𝑙_𝑖=1𝑞𝐙^supplier_calc^→𝑖_𝑚≤∑^𝑙_𝑖=1𝑞^supplier_meas^𝑖_𝑚+Δ,if∑^𝑙_𝑖=1𝑞^supplier_meas^𝑖_𝑚≥−∑^𝑘_𝑗=1𝑞^consumer_meas^𝑗_𝑚,𝐹_scenario=11or21or31_𝑙,_𝑖∑₌₁𝑞^supplier_meas^𝑖_𝑚−Δ≤∑^𝑙_𝑖=1𝑞^supplier_calc^→𝐙^𝑖_𝑚≤−∑^𝑘_𝑗=1𝑞^consumer_meas^𝑗_𝑚+Δ,if∑^𝑙_𝑖=1𝑞^supplier_meas^𝑖_𝑚<−∑^𝑘_𝑗=1𝑞^consumer_meas^𝑗_𝑚,𝐹_scenario=11or21or31_𝑘,_𝑗∑₌₁𝑞𝐙^consumer_calc^→−𝑞^consumer_meas^𝑗_𝑚²−𝜀≤0,if𝐹_scenario=13or23or33,_𝑙𝑗∑₌₁𝑞𝐙^supplier_calc^→−𝑞^supplier_meas^𝑗_𝑚²−𝜀≤0,otherwise. (3.8)(3.9) For numerical solution of the problem (3.5)–(3.9), the following default values for the second and third groups of constraints in (3.7) are recommended: 𝑔_𝑞^𝑠_𝑚=max−𝜀;̃min𝑞^supplier_meas^𝑠_𝑚−𝑄_discrep;𝑞^supplier_meas^𝑠_𝑚−Δ𝑞_𝑚,𝑠=𝑔1,𝑙,_𝑞𝑠𝑚=min𝑞^consumer_meas𝑚^𝑠−𝑙−𝑄_discrep;𝑞^consumer_meas𝑚^𝑠−𝑙−Δ𝑞_𝑚,𝑠=𝑓𝑙+1,𝑙+𝑘,_𝑞𝑠𝑚=max𝑞^supplier_meas^𝑠_𝑚+𝑄_discrep;𝑞^supplier_meas^𝑠_𝑚+Δ𝑞_𝑚,𝑠=𝑓1,𝑙,_𝑞𝑠𝑚=min𝜀;̃max𝑞^consumer_meas𝑚^𝑠−𝑙+𝑄_discrep;𝑞^consumer_meas𝑚^𝑠−𝑙+Δ𝑞_𝑚,𝑠=𝑙+1,𝑙+𝑘,𝑚=0,𝑁_𝑡. (3.10)

In practice, quite a suitable procedure for numerical solution of the problem (3.5)–(3.9) at the 𝑚-th time step is the well-known method of modified Lagrange functions [31–33]. This method provides for constructing a modified Lagrange function defined as 𝐿_𝑐^→𝝁_𝑟,^→𝐙_𝑚=𝑃𝐹_scenario,^→𝐙_𝑚+0.5𝑐_𝑟^{−1𝑙+𝑘+1}_𝑗=1max0;𝜇^𝑗_𝑟+𝑐_𝑟𝑔_𝑗^→𝐙_𝑚²−𝜇^𝑗_𝑟², (3.11) where ^→𝐠(^→𝐙_𝑚)∈𝑅^𝑙+𝑘+1 is a composite function of constraints comprising the second and third groups of constraints in (3.7), and a constraint given by the inequality (3.9). For the Lagrange multiplier vector ^→𝝁_𝑟 given at the 𝑟th iteration and the value of the scalar penalty parameter 𝑐_𝑟, the vector 𝐙^→_𝑚 is defined as a minimum of the function (3.11) with simple constraints on the variables (^→𝐠_𝑚≤^→𝐙_𝑚≤^→𝐟_𝑚) (see the first group of constraints in (3.7)). The problem of minimum search for the function (3.11) with simple constraints on variables can be solved by a modified conjugate direction method [31, 34, 35], which is stable with respect to the accumulation of arithmetic errors. Next, we calculate: 𝜇^𝑗_𝑟+1=max0;𝜇^𝑗_𝑟+𝑐_𝑟𝑔_𝑗^→𝐙^𝑟_𝑚,𝑗=𝑐1,𝑙+𝑘+1,_𝑟+1=𝛽𝑐_𝑟as^→𝐠^→𝐙^𝑟_𝑚0>𝛾^→𝐠^→𝐙_𝑚^𝑟−1₀,𝑐_𝑟,otherwise. (3.12) The initial vector ^→𝝁₀ is chosen as close as possible to the optimal vector ^→𝝁_opt. 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) with simple constraints artificially ill conditioned.

The successive solution of the problem (3.5)–(3.9) at the (𝑁_𝑡+1) time steps makes it possible to recover the agreed Dirichlet boundary conditions for all the pipeline system boundaries within the chosen computational scenarios (see above). Upon their recovery based on the discrete values obtained, it is reasonable to interpolate the boundary conditions. Cubic spline interpolation performs well in this case [28].

In order to obtain spacetime distribution of the recovered transient nonisothermal gas mixture flow parameters, one should conduct a numerical analysis of the model (2.1)–(2.13) closed by interpolated boundary and initial conditions. It should be noted here that in practice, if there is no information on initial conditions, they can be approximated by quasi-steady-state simulation results at the zero time layer.

4. On the Criterion in the Comparative Analysis of Finding Solutions

The above approach to the numerical recovery of gas dynamic parameters of gas mixture flows through pipeline systems based on full-scale measurements gives a number of alternative solutions. This is associated, first of all, with a set of computational scenarios involved and ambiguity of building the vector-function of controlled boundary conditions. In order to choose the best approximation of space-time distributions of real flow parameters, one should propose a criterion to compare the calculated gas dynamic parameters. Such a criterion can be developed by quantitative assessment of the fit between calculated and measured parameters of gas mixture pressure versus time at each IP.

For this purpose, let us introduce the so-called identification factor in the first sense for the 𝑗th IP: Ident_Level_1_𝑗=Δ𝜏⁻¹𝜕𝑝𝜕𝑡^IPcalc(𝑡)−𝜕𝑝𝜕𝑡^IPmeas(𝑡)^𝑗₁=Δ𝜏⁻¹_Δ𝜏𝜕𝑝𝜕𝑡^IPcalc(𝑡)−𝜕𝑝𝜕𝑡^IPmeas(𝑡)_𝑗𝑑𝑡,𝑗=1,𝑀_IP. (4.1) The identification factor in the second sense in our case is written as Ident_Level_2_𝑗=Δ𝜏⁻¹𝑝^IP_calc(𝑡)−𝑝^IP_meas(𝑡)^𝑗₁=Δ𝜏⁻¹_Δ𝜏𝑝^IP_calc(𝑡)−𝑝^IP_meas(𝑡)_𝑗𝑑𝑡,𝑗=1,𝑀_IP. (4.2)

Note that for the time interval Δ𝜏, the level of identification of actual time histories of physical gas flow parameters by calculated time histories in the small neighborhood of 𝑗-th IP, other conditions being equal, will be the higher, the smaller is the value of the corresponding identification factor in the first and second sense.

For simultaneous assessment of identification level in the first and second senses, the law of conventional coloring of IPs is used subject to achieved in its neighborhood identification level. One will take the small neighborhood of 𝑗th IP for the time interval Δ𝜏 be at high identification level in the first and second senses, if at least one of the following two conditions is satisfied: (1)Ident_Level_1_𝑗<𝐶^min_1_BlueIP; (2)Ident_Level_2_𝑗≤𝐶^min_2_GreenIP; (3)Ident_Level_1_𝑗≤𝐶^max_1_BlueIP and 𝐶^min_2_GreenIP<Ident_Level_2_𝑗<𝐶^min_2_BlueIP. In order to make the practical results more demonstrative, the IP of interest will be denoted by a green circle because of the law of coloring usage in the IP layout diagram and the identification level in its neighborhood will be called green identification level.

In a similar situation, the identification level is considered satisfactory in the first and second senses, if at least one of the following two conditions is satisfied: (1)𝐶^max_1_BlueIP<Ident_Level_1_𝑗 and 𝐶^min_2_GreenIP<Ident_Level_2_𝑗<𝐶^min_2_BlueIP; (2)𝐶^min_1_BlueIP≤Ident_Level_1_𝑗≤𝐶^max_1_BlueIP and 𝐶^min_2_BlueIP≤Ident_Level_2_𝑗≤𝐶^max_2_BlueIP. In this case, the IP of interest in the IP layout diagram will be denoted by a blue circle (blue identification level).

Achievement of the identification level disputable in the first and second sense (orange identification level) is characterized by satisfying simultaneously the following two conditions: (1)𝐶^min_1_BlueIP≤Ident_Level_1_𝑗≤𝐶^max_1_BlueIP; (2)𝐶^max_2_BlueIP<Ident_Level_2_𝑗. As a rule, such IPs display a systematic error in gas mixture flow pressure measurements. Such IPs need to be carefully analyzed by specialists operating the simulated pipeline system.

If the above combinations of conditions are not satisfied, a conclusion is drawn that there is no identification in the first and second senses in the small neighborhood of this IP for the time interval Δ𝜏. In this case, the IP will be denoted by a red circle in the IP layout diagram, and the lack of identification in its neighborhood corresponds to the red level.

Analysis of the fit between time histories in the first and second sense does not allow us to account for the influence of individual discrepancy spikes in measurement results on the assessment of the achieved identification level to the full extent. Therefore, additional analysis of the fit between calculated and measured functions in the third sense is required in the neighborhood of the 𝑗th IP as follows: Ident_Level_3_𝑗=0,ifIdent_Level_1_𝑗≤𝐶^max_1_BlueIP;𝑝^IP_calc(𝑡)−𝑝^IP_meas(𝑡)₀=sup_Δ𝜏𝑝^IP_calc(𝑡)−𝑝^IP_meas(𝑡)_𝑗,otherwise. (4.3) According to the method of comparison of calculated gas dynamic parameters described here, the identification level established in the first and second senses should be lowered: from green to blue, if 𝐶^min_3_GreenIP<Ident_Level_3_𝑗≤𝐶^sup_3_IP; from blue to orange, if 𝐶^min_3_BlueIP<Ident_Level_3_𝑗≤𝐶^sup_3_IP; from orange to red, if 𝐶^min_3_OrangeIP<Ident_Level_3_𝑗≤𝐶^sup_3_IP; from any color to red, if 𝐶^sup_3_IP<Ident_Level_3_𝑗. The procedure of lowering the identification level for each IP can be done only once, that is, successive lowering of the green level to the blue one, the blue one, to the orange, and the orange, to the red is not permitted.

The overall assessment of the actual identification level achieved by the 𝑟-th computational gas dynamic mode of actual gas mixture flow through the pipeline system of interest over the given time interval is done using the following formula: 𝑃_{_Ident𝑟}=𝑆_green𝑃_Ident_green𝑟+𝑆_blue𝑃_Ident_blue𝑟+𝑆_orange𝑃_Ident_orange𝑟×𝑆_green𝐶^min_2_GreenIP⁻¹𝑀_IP⁻¹,𝑟=1,𝑉_CFD,𝑃_Ident_green𝑟=0,if𝐿^𝑟_GreenIP=0,_{𝐿𝑟GreenIP𝑗}∑₌₁max𝐶^min_2_GreenIP,Ident_Level_2_𝑗𝑟⁻¹,otherwise,𝑃_Ident_blue𝑟=0,if𝐾^𝑟_BlueIP=0,_{𝐾𝑟BlueIP𝑗}∑₌₁max𝐶^min_2_BlueIP,Ident_Level_2_𝑗𝑟⁻¹,otherwise,𝑃_Ident_orange𝑟=0,if𝑁^𝑟_OrangeIP=0,_{𝑁𝑟OrangeIP𝑗}∑₌₁Ident_Level_2_𝑗𝑟⁻¹,otherwise. (4.4)

The solution to the problem of numerical recovery of gas flows in the simulated pipeline system will be a unique identified gas flow (IGF) that a priori satisfies the defined requirements (constraints) and is compliant with the pipeline system’s real physics of accident-free operation and characterized by the highest value of the quantitative index of identification level (4.4). Thus, the relation 𝑃_Ident_{CFD_ID}=max_{1≤𝑟≤𝑉CFD}𝑃_Ident_𝑟. (4.5) It should be emphasized that the IGF status is assigned to the gas dynamic flow developed only if the following inequality is true for the corresponding prevalence factor of green, blue and orange IPs 𝐹_{CFD_ID}: 𝐹_{CFD_ID}=𝑀⁻¹_IP𝐿^CFD_ID_GreenIP+𝐾^CFD_ID_BlueIP+𝑁^CFD_ID_OrangeIP≥𝐶_{CFD_ID}. (4.6) Otherwise, a conclusion is drawn that the required identification level was not achieved and the IGF was not found, that is, the recovery problem for the gas flows in the gas pipeline system was not solved.

5. On Arrangement of Problem Solution Process

General principles of numerical implementation of identification problem (3.1)–(3.4) can be illustrated by flow diagram for the numerical implementation of the problem (Figure 2). Assignment of new 𝐙^→ and new value of flag 𝐹_scenario is made by users of the method independently proceeding from their research experience of such pipeline systems. As a rule, in practice𝐙^→ is fixed from one time layer to another. This fact is reflected in diagram (see Figure 2). However, it should be noted that possibility of its modification under its layer wise solution was stipulated in problem statement (3.5)–(3.9) (as 𝐙^→_𝑚 is used in (3.5)–(3.9)).

Figure 2 

Flow diagram for the numerical implementation of the problem.

Unique IGF selection from several versions of calculated gas dynamic parameters filed to database of calculation results (DCR) allows compensating negative influence of random measurement errors on credibility of numerical recovery of transport flows.

6. Results of Practical Application

Efficiency of the method of numerical recovery of gas flows in trunkline systems proposed in the paper was demonstrated in 2010–2012 in production simulations at GAZPROM Mezhregiongaz Moscow LLC within the Alfargus/Mosregiongaz Computer System (Figure 3).

Figure 3 

Example of the Alfargus/Mosregiongaz computer system application in the control room of GAZPROM Mezhregiongaz Moscow LLC.

The method was used for numerical recovery of the flow of natural gas delivered (from a single supplier) to consumers through seven branches of the Moscow Gas Ring (MGR). MGR has a total length of over 200 km and more than 130 consumer branches. The flow was recovered at 106 IPs, which were relatively uniformly distributed over the gas pipeline ring. In Figure 4, the given IPs were indicated in the form of circles of different diameter and color depending on its arrangement on the boundaries of “supplier branches” (dark red circles), and on the boundaries of “consumer branches” (pink circles) and in network and on the boundaries of uncontrolled by flow meters branches (yellow circles).

Figure 4 

Example of IP distribution in the MRG model.

The transport flow is transient nonisothermal gas flow. The example of flow diagram (i.e., recovered flow direction and numerical estimates of volumetric flow rate of natural gas (dimension: thousand cubic meters per day) in accordance with color gradation) in the South-East MRG sector (temporal section) was shown in Figure 4. In table on the right of Figure 5 one can see quantitative estimates of gas flow rate distribution (column 2, dimension: thousand standard cubic meters per day) and gas pressure (column 3, dimension: gauge atmospheres) for recovered flows in specific branches in the South-East MRG sector (temporal section). In the first column of the table under consideration description of branches are given in topographical map reference. The example of diagram correlation of time history of calculated and measured estimates of pressure and mass flow rates for one from the IPs, which is used in MGR (gas flow temperature was measured with a poor accuracy and long time intervals and was not suitable for comparative analysis) was shown in Figure 6. It should be noted that measurement results underwent preliminary verification and smoothing. The recovered gas flow parameters were used to analyze the performance of MGR, and to detect and localize the sources of discrepancy in estimated volumes of gas supply through MGR [36].

Figure 5 

Example of flow diagram in the South-East MRG sector (temporal section).

(a)

(b)

Figure 6 

Example of curve correlation of calculated and measured pressure history (a) and mass flow rate (b) for one from the IPs used in MGR.

As a demonstration of natural gas pressure, volumetric gas flow rate and temperature distribution along the length of MGR circular segment (see Figure 4) for one time layer were given on Figure 7. Fragmentary nature of flow rates space distribution containing breaks made for gas extraction and discharge in MGR. In continuous space distribution of transported gas pressure the extremes basically made for gas discharge to MGR circular segment by different suppliers and gas extraction by major consumers. Certain breaks in temperature distribution along the length of MGR circular segment are explained by branches existence and multiline structure of MGR circular segment—for some segments different pipeline threads are shown.

(a)

(b)

Figure 7 

Example of space distributions “pressure-volumetric flow rate” (a) and “pressure-temperature” (b) along the length of MGR circular segment (temporal section).

Earlier versions of the flow recovery method were used to investigate trunkline accidents and to train gas pipeline operators in efficient pipeline control under conditions as close as possible to real operation of gas transmission and delivery systems using high-accuracy computer simulators [16, 37].

7. Conclusion

The paper describes a version of the practical method for numerical recovery of transient gas flow in gas trunkline systems by setting up and solving a special identification problem. Simulations of gas mixture flow parameters in extended branched pipelines, involving calculations of the target function and constraint function for the identification problem of interest, are done in the 1D statement. The method can be implemented on computers available to a wide range of pipeline industry specialists. Efficiency of this method as applied to practical simulations was demonstrated, when it was used to provide computerized support of operator decisions in gas distribution companies, investigate accidents in gas pipeline systems and train gas pipeline operators in pipeline control under normal and contingency conditions using research and high-accuracy computer simulators.

Nomenclature

𝐶_{CFD_ID}: Empirical constant, the value of which is chosen based on the experience of doing simulations with the recovery method described here, 𝐶_{CFD_ID}>0 CFD_ID: Index, will be true for the identified gas flow𝐶^min_2_GreenIP, 𝐶^min_1_BlueIP, 𝐶^max_1_BlueIP, 𝐶^min_2_BlueIP, 𝐶^max_2_BlueIP, 𝐶^min_3_GreenIP, 𝐶^min_3_BlueIP, 𝐶^min_3_OrangeIP, 𝐶^sup_3_IP: Given empirical constants 𝑐_𝑟: Scalar penalty parameter 𝐷_𝑚: Binary diffusion coefficient of the 𝑚th component 𝑓: Gas pipeline flow section area 𝐹_{CFD_ID}: Corresponding prevalence factor of green, blue and orange IPs 𝐹_scenario: Flag of involved computational scenarios of the identification problem 𝑔: Gravity acceleration modulus ^→𝐠(𝑡)∈𝑅^𝑛,^→𝐟(𝑡)∈𝑅^𝑛: Given vector-functions that establish limits in simple constraints on the vector-function of controlled boundary conditions based on structural and operational features of the simulated pipeline system, ^→𝐠(𝑡)<^→𝐟(𝑡) ^→𝐠_𝑞(𝑡)∈𝑅^𝑙+𝑘,^→𝐟_𝑞(𝑡)∈𝑅^𝑙+𝑘: Given vector-functions that establish limits in constraints providing a-priory preservation of the defined purpose of the branch over the time interval Δ𝜏 (i.e., a gas supplier cannot become a gas consumer and vice versa, see above), ^→𝐠_𝑞(𝑡)<^→𝐟_𝑞(𝑡) ^→𝐠(^→𝐙_𝑚)∈𝑅^𝑙+𝑘+1: Composite function Ident_Level_1_𝑗: So-called identification factor in the first sense for the 𝑗th Identification Point (IP) Ident_Level_2_𝑗: So-called identification factor in the second sense for the 𝑗th IP Ident_Level_3_𝑗: So-called identification factor in the third sense is required in the neighborhood of the 𝑗th IP ^→(𝑛)𝐢: Unit vector of the 𝑂𝑥 axis of the 𝑛th pipeline 𝐾^𝑟_BlueIP: Number of blue IPs for the 𝑟th computational gas dynamic mode 𝐾: Heat conductivity of the gas mixture 𝑘: Given number of consumer branches  𝐿: Identifier of the considered pipeline cross sections 𝐿₀,𝐿₁,𝐿₂: Metrics 𝐿_𝑐(^→𝝁_𝑟,𝐙^→_𝑚): Modified Lagrange function 𝐿^𝑟_GreenIP: Number of green IPs for each 𝑟th computational gas dynamic mode 𝑙: Given number of supplier branches, through which the gas mixture enters the pipeline system over the time Δ𝜏 𝑀: Identifier of the considered pipeline cross sections 𝑀_difference: Number of analyzable pressure differences between neighbor IPs 𝑀_IP: Number of IPs 𝑁: Number of pipes constituting the simulated joint 𝑁_𝑡: Number of time layers separated from each other by a given uniform step Δ𝜏 𝑁_𝑆: Number of mixture components 𝑁^𝑟_OrangeIP: Number of orange IPs for the 𝑟th computational gas dynamic mode ^→(0)𝐧: Unit vector of the normal to the volume ⁽⁰⁾𝑉 𝑃[𝑡,𝐹_scenario,^→𝐙(𝑡)]: Target function 𝑃_Ident_𝑟: Actual identification level achieved by the 𝑟th computational gas dynamic mode of actual gas mixture flow through the pipeline system 𝑃_Ident_{CFD_ID}: Quantitative index of identification level 𝑝: Static pressure in the gas mixture 𝑝_Joint: Static pressure of the gas mixture in the joint (i.e., in the inner space of the volume ⁽⁰⁾𝑉 ^→𝐩^IP_calc[𝑡,^→𝐙(𝑡)]∈𝑅^𝑀IP: Vector-function describing the time variation of calculated estimates of gas mixture pressure at the IP ^→𝐩^IP_meas(𝑡)∈𝑅^𝑀IP: Given vector-function describing the time variation of measured estimates of gas mixture pressure at the IP 𝑄: Specific (per unit volume) power of heat sources 𝑄_discrep: Predefined empirical constant corresponding to the minimum value of the module of total discrepancy between estimated gas supply volumes, which is of practical significance for the analysis of the simulated gas pipeline network ^→𝐪^consumer_calc[𝑡,𝐙^→(𝑡)]∈𝑅^𝑘: Vector-function that describes the time variation of calculated estimates of gas mixture mass flow rate at outlets of consumer branches ^→𝐪^supplier_calc[𝑡,^→𝐙(𝑡)]∈𝑅^𝑙: Vector-function describing the time variation of calculated estimates of mass flow rates of the gas mixture at outlets of supplier branches ^→𝐪^consumer_meas(𝑡)∈𝑅^𝑘: Given vector-function that describes the time variation of measured or declared estimates of gas mass flow rates at outlets of consumer branches ^→𝐪^supplier_meas(𝑡)∈𝑅^𝑙: Given vector-function describing the time variation of measured estimates of mass flow rates of the gas mixture at outlets of supplier branches 𝑅=𝑓/𝜋: Inner pipe radius 𝑅^𝑛: n-dimensional Euclidean space 𝑟: Iteration index of the method of modified Lagrange functions 𝑆_green, 𝑆_blue, 𝑆_orange: Scalar weight factors used when establishing quantitative indices of the identification level achieved, (𝑆_green>𝑆_blue≫𝑆_orange>0) 𝑠: Number of outlets of consumer branches, 𝑠≤𝑘 ^(𝑛)𝑠: Auxiliary functions  𝑇: Gas mixture temperature 𝑇_am: Space-time distribution of ambient temperature at the domain boundary 𝑡: Time (marching variable) ⁽⁰⁾𝑉: Vanishingly small volume of the simulated joint of trunklines;  𝑉_CFD: number of alternative computational gas dynamic modes obtained by practical implementation of Section 4 provisions 𝑤: Projection of the vector of gas velocity averaged over the pipeline cross section to the pipeline’s geometric axis of symmetry (on the assumption of developed flow turbulence) 𝑥: Spatial coordinate along the geometric pipeline axis (spatial variable) 𝑌_𝑚=𝜌_𝑚/𝜌: Relative mass fraction of the 𝑚th component 𝑌_𝑚,Joint: Relative mass fraction of the 𝑚th component in the joint 𝑦_𝑖: ith component of 𝐘^→∈𝑅^𝑛 𝐙^→(𝑡): Vector-function of controlled boundary conditions that describes variations in boundary conditions at pre-selected pipeline boundaries during calculation of the identification problem 𝐙^→: Vector of controlled boundary conditions that describes variations in boundary conditions at pre-selected pipeline boundaries during calculation of the identification problem 𝑧₁: Spatial coordinate of the point in the pipeline axis reckoned from an arbitrary horizontal plane vertically upward (for natural gas trunklines, along the Earth radius) 𝛽: Given integer-valued parameter multiple of 10 ^(𝑛)𝛾: Parameter, 𝑛=1,𝑁 Δ=𝑄_discrep/2 Δ𝑞(𝑡): Given function that describes the time variation of the module of declared discrepancy between estimated gas supply volumes Δ𝑡: Preceding time interval ^(𝑛)Δ𝑋: Elementary segment between the cross sections ^(𝑛)𝑓_𝐿 and ^(𝑛)𝑓_𝑀 bounding each elementary volume ^(𝑛)𝑉,𝑛=1,𝑁 Δ𝑋: Comparison parameter, Δ𝑋→0 Δ𝜏: Given time interval 𝜀: Specific (per unit mass) internal energy of the mixture 𝜀_𝑚: Specific (per unit mass) internal energy of the 𝑚th component 𝜀_Joint: Specific (per unit mass) internal energy is of the gas mixture in the joint 𝜀: Predefined small quantity that establishes the minimum difference between calculated and measured mass flow rates in the second sense of close fit, 𝜀>0 𝜀̃: Predefined small quantity to establish a boundary of practical significance of the gas mixture mass flow rate through the branch (a lower-modulus gas mass flow rate is considered zero), 𝜀̃>0 ^(𝑛)Θ: Auxiliary function 𝜆: Hydraulic friction coefficient in the widely known Darcy-Weisbach formula ^→𝝁_𝑟∈𝑅^𝑙+𝑘+1: Vector of Lagrange multipliers ^→𝝁₀: Initial vector of Lagrange multipliers ^→𝝁_opt: Optimal vector of Lagrange multipliers 𝚷(𝑡): Constraint set 𝜋: Pythagorean number 𝜌: Gas mixture density 𝜌_Joint: Gas mixture density in the joint 𝜌_𝑚: Reduced density of the 𝑚th component (mass of the 𝑚th component in unit volume of the mixture) Φ(𝑇,𝑇_am): Function characterizes the heat exchange of the gas flow core through the boundary gas layer, pipe wall and insulation with the environment 𝜒: Perimeter of the cross section having an area of 𝑓 Ω_𝑚(𝐹_scenario): Auxiliary multiplier 𝝎^→_I-II[𝑡,^→𝐙(𝑡)]∈𝑅^{𝑀difference}: Auxiliary vector-function; 𝑀_difference is the number of given pairs of IPs that determine the controlled natural and virtual pressure drops in the gas mixture ⟨unequality(Δ𝜏,𝐹_scenario,^→𝐙(𝑡))⟩: Formal representation of an additional limiting inequality {𝑆_mix}: Formal expression corresponding to a set of parameters defining the described quantity (^→(0)𝐧⋅^→(𝑛)𝐢): Scalar product of the vectors ^→(0)𝐧and ^→(𝑛)𝐢 Left superscript “(𝐧)” denotes that some parameter belongs to pipe number 𝑛.