# Numerical Recovery of Gas Flows in Pipeline Systems

Authors: Vadim E. Seleznev
Journal: Journal of Applied Mathematics
 (2012)

---

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


---

## Body

## 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 ofN pipelines.(i) for each pipe:(2.1)∂(ρf)∂t+∂(ρwf)∂x=0,(2.2)∂(ρYmf)∂t+∂(ρYmwf)∂x-∂∂x(ρfDm∂Ym∂x)=0,m=1,NS-1¯,YNS=1-∑m=1NS-1Ym,(2.3)∂(ρwf)∂t+∂(ρw2f)∂x=-f(∂p∂x+gρ∂z1∂x)-π4λρw|w|R,(2.4)∂∂t[ρf(ε+w22)]+∂∂x[ρwf(ε+w22)]=-∂(pwf)∂x-ρwfg∂z1∂x-p∂f∂t+Qf+∂∂x(Kf∂T∂x)-Φ(T,Tam)+∂∂x(ρf∑m=1NS{εmDm∂Ym∂x});(ii) junction conditions for each joint:(2.5)∑n=1N([wf]Ls)(n)=0;∑n=1N([f∂T∂x]Ls)(n)=0,∑n=1N([f∂Ym∂x]Ls)(n)=0,m=1,NS-1¯,(2.6)TL(n)=TL(ξ),εJoint=εL(n),(εm)L(n)=(εm)L(ξ),ρJoint=ρL(n),pJoint=pL(n),kL(n)=kL(ξ),(Dm)L(n)=(Dm)L(ξ),Ym,Joint=(Ym)L(n),(z1)L(n)=(z1)L(ξ)∀n,ξ∈1,N¯and∀m∈1,NS¯;(iii) fitting conditions for each joint:(2.7)∂ρJoint∂t+∑n=1NΘ(n)[∂(ρw)∂x]L(n)=0,(2.8)∂(ρYm)Joint∂t+∑n=1NΘ(n)[∂(ρYmw)∂x-∂∂x(ρDm∂Ym∂x)]L(n)=0,YNS,Joint=1-∑m=1NS-1Ym,Joint,m=1,NS-1¯,(2.9)[∂(ρw)∂t+∂(ρw2)∂x+∂p∂x+gρ∂z1∂x+0.25λρw|w|R-1]L(n)=0,n=1,N¯,(2.10)∂(ρε)Joint∂t+∑n=1NΘ(n)[∂(ρεw)∂x+p∂w∂x-0.25λρ|w|3R-1-∂∂x(K∂T∂x)]L(n)+∑n=1NΘ(n)[Φ(T,Tam)f-∑m=1NS∂∂x(ρεmDm∂Ym∂x)-Q]L(n)=0;(iv) EOS and additional relations:(2.11)p=p({Smix}),ε=ε({Smix}),K=K({Smix}),T1=T2=⋯=TNS=T,εm=εm({Smix}),Dm=Dm({Smix}),m=1,NS¯;(v) auxiliary geometric relations:(2.12)s(n)=-(n→(0)·i→(n))={1,if(n→(0)·i→(n))<0,-1,if(n→(0)·i→(n))>0,(2.13)Θ(n)=V(n)V=(γf)(n)∑k=1N(γf)(k),0<Θ(n)<1,∑l=1NΘ(l)=1;γ(n)=ΔX(n)ΔX-,n=1,N¯;ΔX(n)=|xL(n)-xM(n)|=s(n)(xL(n)-xM(n))=[s(xL-xM)](n).
The function Φ(T,Tam) 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 f from the transported gas to the environment (Φ(T,Tam)>0 corresponds to heat removal; Tam 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 γ,n(n)=1,N¯ (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 γ,n(n)=1,N¯, which makes it ambiguous. However, one can demonstrate that this model remains correct for arbitrary values ofγ,n(n)=1,N¯. 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 L1 (in our case, it is defined by means of the octahedral vector norm, that is, L1=∥Y→∥1=∑i=1n|yi|,Y→∈Rn) or metric L2 (in our case, it is defined by means of the Euclidean vector norm,that is, L2=∥Y→∥2=[∑i=1nyi2]0.5,Y→∈Rn); (3) close fit between two functional relations within their uniform deviation, that is, in the metric L0 (in our case, it is defined by means of the cubic vector norm, that is, L0=∥Y→∥0=max1≤i≤n|yi|,Y→∈Rn).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:(3.1)P[t,Fscenario,Z→(t)]→minZ→(t)∈Π(t)⊂Rn.In our case, components of the vector-functionZ→(t) 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 u(u≤l) outlets of supplier branches and at s(s≤k) outlets of consumer branches, where k 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 (l-u) outlets of supplier branches and at (k-s) outlets of consumer branches. Hence, the total number of variable components in the boundary conditions is n=l+k. 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, n<l+k. 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 problemFscenario 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 Fscenario will be demonstrated below (see (3.2)–(3.9)).The target functionP[t,Fscenario,Z→(t)] 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:(3.2)P[t,Fscenario,Z→(t)]={∥p→calcIP[t,Z→(t)]-p→measIP(t)∥2,ifFscenario<21,∥∂∂tp→calcIP[t,Z→(t)]-∂∂tp→measIP(t)∥2,ifFscenario<31,∥ω→I-II[t,Z→(t)]∥2,otherwise,ω→I-II[t,Z→(t)]={p→calcIP[t,Z→(t)]-p→measIP(t)}I+{p→measIP(t)-p→calcIP[t,Z→(t)]}II.
The vector-function p→calcIP[t,Z→(t)] 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 Z→(t).The first form of the target function (i.e., forFscenario<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 (forFscenario<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 Mdifference≫MIP 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):(3.3)Z→(t)∈Π(t)={Z→(t)∈Rn:g→(t)≤Z→(t)≤f→(t),[gq(t)]s≤[qcalcsupplier[t,Z→(t)]]s≤[fq(t)]s,s=1,l¯,[gq(t)]s≤[qcalcconsumer[t,Z→(t)]]s-l≤[fq(t)]s,s=l+1,l+k¯,〈unequality(Δτ,Fscenario,Z→(t))〉}.
The vector-function q→calcconsumer[t,Z→(t)] 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 Z→(t). The following law of signs is true for the vector-function q→calcconsumer[t,Z→(t)]: if [qcalcconsumer[t,Z→(t)]]i<0, the gas mixture moves from the gas pipeline system to the consumer, i=1,k¯. 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:
(3.4)〈unequality(Δτ,Fscenario,Z→(t))〉={-∫Δτ∑j=1k[qmeasconsumer(t)]jdt-Δ≤∫Δτ∑i=1l[qcalcsupplier[t,Z→(t)]]idt≤∫Δτ∑i=1l[qmeassupplier(t)]idt+Δ,if∫Δτ∑i=1l[qmeassupplier(t)]idt≥-∫Δτ∑j=1k[qmeasconsumer(t)]jdt,(Fscenario=11or21or31),∫Δτ∑i=1l[qmeassupplier(t)]idt-Δ≤∫Δτ∑i=1l[qcalcsupplier[t,Z→(t)]]idt≤-∫Δτ∑j=1k[qmeasconsumer(t)]jdt+Δ,if∫Δτ∑i=1l[qmeassupplier(t)]idt<-∫Δτ∑j=1k[qmeasconsumer(t)]jdt,(Fscenario=11or21or31),∥q→calcconsumer[t,Z→(t)]-q→measconsumer(t)∥2-ε≤0,ifFscenario=13or23or33,∥q→calcsupplier[t,Z→(t)]-q→meassupplier(t)∥2-ε≤0,otherwise.
The vector-function q→calcsupplier[t,Z→(t)] 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. [qcalcsupplier[t,Z→(t)]]i>0 or [qmeassupplier(t)]i>0 means that the gas mixture enters the gas pipeline system, i=1,l¯. The following law of signs is true for components of the vector-function q→measconsumer(t): if [qmeasconsumer(t)]i<0, the gas mixture moves from the gas pipeline system to the consumer, i=1,k¯.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 (Nt+1) time layers separated from each other by a given uniform step Δτ. The m=0 layer will correspond to the lower boundary of the time interval Δτ, and the m=Nt 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 Δt. 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:
(3.5)P(Fscenario,Z→m)→minZ→m∈Πm⊂Rn;(3.6)P(Fscenario,Z→m)={∑i=0MIP([pcalcIP(Z→)-pmeasIP]mi-Ωm(Fscenario)[pcalcIP(Z→)-pmeasIP]m-1i)2,ifFscenario<31;∑i=0Mdifference([pcalcIP(Z→)-pmeasIP]mI,i+[pmeasIP-pcalcIP(Z→)]mII,i)2,otherwise,(3.7)Z→m∈Πm={Z→m∈Rn:g→m≤Z→m≤f→m,[gq]ms≤[qcalcsupplier(Z→)]ms≤[fq]ms,s=1,l¯,[gq]ms≤[qcalcconsumer(Z→)]ms-l≤[fq]ms,s=l+1,l+k¯,〈unequality(Fscenario,Z→)〉m},m=0,Nt¯,
where
(3.8)Ωm(Fscenario)={0,ifm=0orFscenario<21,1,otherwise,(3.9)〈unequality(Fscenario,Z→)〉m={-∑j=1k[qmeasconsumer]mj-Δ≤∑i=1l[qcalcsupplier(Z→)]mi≤∑i=1l[qmeassupplier]mi+Δ,if∑i=1l[qmeassupplier]mi≥-∑j=1k[qmeasconsumer]mj,(Fscenario=11or21or31),∑i=1l[qmeassupplier]mi-Δ≤∑i=1l[qcalcsupplier(Z→)]mi≤-∑j=1k[qmeasconsumer]mj+Δ,if∑i=1l[qmeassupplier]mi<-∑j=1k[qmeasconsumer]mj,(Fscenario=11or21or31),∑j=1k([qcalcconsumer(Z→)-qmeasconsumer]mj)2-ε-≤0,ifFscenario=13or23or33,∑j=1l([qcalcsupplier(Z→)-qmeassupplier]mj)2-ε-≤0,otherwise.
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:
(3.10)[gq]ms=max{-ε~;min([qmeassupplier]ms-Qdiscrep;[qmeassupplier]ms-Δqm)},s=1,l¯,[gq]ms=min{[qmeasconsumer]ms-l-Qdiscrep;[qmeasconsumer]ms-l-Δqm},s=l+1,l+k¯,[fq]ms=max{[qmeassupplier]ms+Qdiscrep;[qmeassupplier]ms+Δqm},s=1,l¯,[fq]ms=min{ε~;max([qmeasconsumer]ms-l+Qdiscrep;[qmeasconsumer]ms-l+Δqm)},s=l+1,l+k¯,[fq]ms=min{ε~;max([qmeasconsumer]ms-l+Qdiscrep;[qmeasconsumer]ms-l+Δqm)},s=m=0,Nt¯.In practice, quite a suitable procedure for numerical solution of the problem (3.5)–(3.9) at the m-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
(3.11)L~c(μ→r,Z→m)=P(Fscenario,Z→m)+0.5cr-1∑j=1l+k+1{[max{0;μrj+crg~j(Z→m)}]2-(μrj)2},
where g~→(Z→m)∈Rl+k+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 μ→r given at the rth iteration and the value of the scalar penalty parameter cr, the vector Z→m is defined as a minimum of the function (3.11) with simple constraints on the variables (g→m≤Z→m≤f→m) (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:
(3.12)μr+1j=max{0;μrj+crg~j(Z→mr)},j=1,l+k+1¯,cr+1={βcras∥g~→(Z→mr)∥0>γ~∥g~→(Z→mr-1)∥0,cr,otherwise.
The initial vector μ→0 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 c0 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 (Nt+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 thejth IP:
(4.1)Ident_Level_1j=Δτ-1∥∂∂tpcalcIP(t)-∂∂tpmeasIP(t)∥1j=Δτ-1∫Δτ|∂∂tpcalcIP(t)-∂∂tpmeasIP(t)|jdt,j=1,MIP¯.
The identification factor in the second sense in our case is written as
(4.2)Ident_Level_2j=Δτ-1∥pcalcIP(t)-pmeasIP(t)∥1j=Δτ-1∫Δτ|pcalcIP(t)-pmeasIP(t)|jdt,j=1,MIP¯.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 j-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 ofjth 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_1j<CBlue IPmin_1; (2)Ident_Level_2j≤CGreen IPmin_2; (3)Ident_Level_1j≤CBlue IPmax_1 and CGreen IPmin_2<Ident_Level_2j<CBlue IPmin_2. 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)CBlue IPmax_1<Ident_Level_1j and CGreen IPmin_2<Ident_Level_2j<CBlue IPmin_2; (2)CBlue IPmin_1≤Ident_Level_1j≤CBlue IPmax_1 and CBlue IPmin_2≤Ident_Level_2j≤CBlue IPmax_2. 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)CBlueIPmin_1≤Ident_Level_1j≤CBlueIPmax_1; (2)CBlueIPmax_2<Ident_Level_2j. 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 thejth IP as follows:
(4.3)Ident_Level_3j={0,ifIdent_Level_1j≤CBlueIPmax_1;∥pcalcIP(t)-pmeasIP(t)∥0=supΔτ|pcalcIP(t)-pmeasIP(t)|j,otherwise.
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 CGreenIPmin_3<Ident_Level_3j≤CIPsup_3; from blue to orange, if CBlueIPmin_3<Ident_Level_3j≤CIPsup_3; from orange to red, if COrangeIPmin_3<Ident_Level_3j≤CIPsup_3; from any color to red, if CIPsup_3<Ident_Level_3j. 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 ther-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:
(4.4)P_Identr={Sgreen[P_Identgreen]r+Sblue[P_Identblue]r+Sorange[P_Identorange]r}×(Sgreen[CGreenIPmin_2]-1MIP)-1,r=1,VCFD¯,[P_Identgreen]r={0,ifLGreenIPr=0,∑j=1LGreenIPr[max{CGreenIPmin_2,Ident_Level_2j}]r-1,otherwise,[P_Identblue]r={0,ifKBlueIPr=0,∑j=1KBlueIPr[max{CBlueIPmin_2,Ident_Level_2j}]r-1,otherwise,[P_Identorange]r={0,ifNOrangeIPr=0,∑j=1NOrangeIPr[Ident_Level_2j]r-1,otherwise.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
(4.5)P_IdentCFD_ID=max1≤r≤VCFD{P_Identr}.
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 FCFD_ID:
(4.6)FCFD_ID=MIP-1(LGreenIPCFD_ID+KBlueIPCFD_ID+NOrangeIPCFD_ID)≥CCFD_ID.
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 Z→ and new value of flag Fscenario is made by users of the method independently proceeding from their research experience of such pipeline systems. As a rule, in practiceZ→ 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 Z→m 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 (Figure3).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 Figure4, 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 Figure4. 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).Example of curve correlation of calculated and measured pressure history (a) and mass flow rate (b) for one from the IPs used in MGR.
(a)
(b)As a demonstration of natural gas pressure, volumetric gas flow rate and temperature distribution along the length of MGR circular segment (see Figure4) 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.Example of space distributions “pressure-volumetric flow rate” (a) and “pressure-temperature” (b) along the length of MGR circular segment (temporal section).
(a)
(b)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

C
CFD
_
ID
: Empirical constant, the value of which is chosen based on the experience of doing simulations with the recovery method described here,CCFD_ID>0
CFD
_
ID: Index, will be true for the identified gas flow
C
Green
IP
min
_
2
,CBlueIPmin_1, CBlueIPmax_1, CBlueIPmin_2, CBlueIPmax_2, CGreenIPmin_3, CBlueIPmin_3, COrangeIPmin_3, CIPsup_3: Given empirical constants
c
r: Scalar penalty parameter
D
m: Binary diffusion coefficient of the mth component
f: Gas pipeline flow section area
F
CFD
_
ID: Corresponding prevalence factor of green, blue and orange IPs
F
scenario: Flag of involved computational scenarios of the identification problem
g: Gravity acceleration modulus
g
→
(
t
)
∈
R
n
,
f
→
(
t
)
∈
R
n: 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, g→(t)<f→(t)
g
→
q
(
t
)
∈
R
l
+
k
,
f
→
q
(
t
)
∈
R
l
+
k: 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), g→q(t)<f→q(t)
g
~
→
(
Z
→
m
)
∈
R
l
+
k
+
1: Composite function
Ident
_
Level
_
1
j: So-called identification factor in the first sense for the jth Identification Point (IP)
Ident
_
Level
_
2
j: So-called identification factor in the second sense for the jth IP
Ident
_
Level
_
3
j: So-called identification factor in the third sense is required in the neighborhood of the jth IP
i
→
(
n
): Unit vector of the Ox axis of the nth pipeline
K
Blue
IP
r: Number of blue IPs for the rth computational gas dynamic mode
K: Heat conductivity of the gas mixture
k: Given number of consumer branches
L: Identifier of the considered pipeline cross sections
L
0
,
L
1
,
L
2: Metrics
L
~
c
(
μ
→
r
,
Z
→
m
): Modified Lagrange function
L
Green
IP
r: Number of green IPs for each rth computational gas dynamic mode
l: Given number of supplier branches, through which the gas mixture enters the pipeline system over the time Δτ
M: Identifier of the considered pipeline cross sections
M
difference: Number of analyzable pressure differences between neighbor IPs
M
IP: Number of IPs
N: Number of pipes constituting the simulated joint
N
t: Number of time layers separated from each other by a given uniform step Δτ
N
S: Number of mixture components
N
Orange
IP
r: Number of orange IPs for the rth computational gas dynamic mode
n
→
(
0
): Unit vector of the normal to the volume V(0)
P
[
t
,
F
scenario
,
Z
→
(
t
)
]: Target function
P
_
Iden
t
r: Actual identification level achieved by the rth computational gas dynamic mode of actual gas mixture flow through the pipeline system
P
_
Iden
t
CFD
_
ID: Quantitative index of identification level
p: Static pressure in the gas mixture
p
Joint: Static pressure of the gas mixture in the joint (i.e., in the inner space of the volume V(0)
p
→
calc
IP
[
t
,
Z
→
(
t
)
]
∈
R
M
IP: Vector-function describing the time variation of calculated estimates of gas mixture pressure at the IP
p
→
meas
IP
(
t
)
∈
R
M
IP: Given vector-function describing the time variation of measured estimates of gas mixture pressure at the IP
Q: Specific (per unit volume) power of heat sources
Q
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
q
→
calc
consumer
[
t
,
Z
→
(
t
)
]
∈
R
k: Vector-function that describes the time variation of calculated estimates of gas mixture mass flow rate at outlets of consumer branches
q
→
calc
supplier
[
t
,
Z
→
(
t
)
]
∈
R
l: Vector-function describing the time variation of calculated estimates of mass flow rates of the gas mixture at outlets of supplier branches
q
→
meas
consumer
(
t
)
∈
R
k: Given vector-function that describes the time variation of measured or declared estimates of gas mass flow rates at outlets of consumer branches
q
→
meas
supplier
(
t
)
∈
R
l: Given vector-function  describing the time variation of measured estimates of mass flow rates of the gas mixture at outlets of supplier branches
R
=
f
/
π: Inner pipe radius
R
n: n-dimensional Euclidean space
r: Iteration index of the method of modified Lagrange functions
S
green, Sblue, Sorange: Scalar weight factors used when establishing quantitative indices of the identification level achieved, (Sgreen>Sblue≫Sorange>0)
s: Number of outlets of consumer branches, s≤k
s
(
n
): Auxiliary functions
T: Gas mixture temperature
T
am: Space-time distribution of ambient temperature at the domain boundary
t: Time (marching variable)
V
(
0
): Vanishingly small volume of the simulated joint of trunklines;
V
CFD: number of alternative computational gas dynamic modes obtained by practical implementation of Section 4 provisions
w: 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)
x: Spatial coordinate along the geometric pipeline axis (spatial variable)
Y
m
=
ρ
m
/
ρ: Relative mass fraction of the mth component
Y
m
,
Joint: Relative mass fraction of the mth component in the joint
y
i: ith component of  Y→∈Rn
Z
→
(
t
): Vector-function of controlled boundary conditions that describes variations in boundary conditions at pre-selected pipeline boundaries during calculation of the identification problem
Z
→: Vector of controlled boundary conditions that describes variations in boundary conditions at pre-selected pipeline boundaries during calculation of the identification problem
z
1: 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
γ
(
n
): Parameter, n=1,N¯
Δ
=
Q
discrep
/
2
Δ
q
(
t
): Given function that describes the time variation of the module of declared discrepancy between estimated gas supply volumes
Δ
t: Preceding time interval
Δ
X
(
n
): Elementary segment between the cross sections fL(n) and fM(n) bounding each elementary volume V(n),n=1,N¯
Δ
X
-: Comparison parameter, ΔX-→0
Δ
τ: Given time interval
ε: Specific (per unit mass) internal energy of the mixture
ε
m: Specific (per unit mass) internal energy of the mth 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
Θ
(
n
): Auxiliary function
λ: Hydraulic friction coefficient in the widely known Darcy-Weisbach formula
μ
→
r
∈
R
l
+
k
+
1: Vector of Lagrange multipliers
μ
→
0: Initial vector of Lagrange multipliers
μ
→
opt: Optimal vector of Lagrange multipliers
Π
(
t
): Constraint set
π: Pythagorean number
ρ: Gas mixture density
ρ
Joint: Gas mixture density in the joint
ρ
m: Reduced density of the mth component (mass of the mth component in unit volume of the mixture)
Φ
(
T
,
T
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 f
Ω
m
(
F
scenario
): Auxiliary multiplier
ω
→
I-II
[
t
,
Z
→
(
t
)
]
∈
R
M
difference: Auxiliary vector-function; Mdifference is the number of given pairs of IPs that determine the controlled natural and virtual pressure drops in the gas mixture
〈
unequality
(
Δ
τ
,
F
scenario
,
Z
→
(
t
)
)
〉: Formal representation of an additional limiting inequality
{
S
mix
}: Formal expression corresponding to a set of parameters defining the described quantity
(
n
→
(
0
)
·
i
→
(
n
)
): Scalar product of the vectors n→(0)and i→(n)
Left superscript“(n)” denotes that some parameter belongs to pipe number n.



---

*Source: 180645-2012-10-02.xml*