# Numerical Adaptation of Pipeline Network Models on Measurement Archive

Authors: Vadim E. Seleznev
Journal: ISRN Applied Mathematics
 (2014)

---

## Abstract

We propose an adaptation method for gas dynamic pipeline network models to enable credible representation of actual properties of real simulation objects. The presentation is illustrated by fitting equivalent pipeline section roughnesses used in the models to accommodate the influence of flow resistance on gas transport parameters. The method is based on the setting up and solution of a series of special parametric identification problems based on a limited set of field measurement data at local (in space) network points. This method can be used by specialists in mathematical modeling of gas transport systems to solve practical parametric identification problems.


---

## Body

## 1. Problem Statement

Under present-day conditions of progress in computational technology, various problems of monitoring and operation of complex pipeline systems are closely related to the problem of numerical recovery of flow parameters in pipelines [1]. Mathematical modeling of fluid transport processes in a pipeline network offers a promising approach to this problem. Various properties of the engineering system of interest are represented by corresponding model parameters. Some properties of a branched pipeline network that are essential for flow parameters may vary with time over its life cycle. If the model parameters corresponding to such properties are static, or if their variation law cannot be defined a priori, then such model parameters should be periodically adjusted or revised. Without loss of generality, this method of pipeline network model adjustment will be presented for the case of numerical fitting of effective (equivalent) roughnesses of inner pipe walls in linear segments of the pipeline system of interest.As we know [2], total losses of pressure Δpfull in viscous matter flowing through pipes to surmount the flow resistance forces are usually defined as a sum of friction losses Δpfric and local losses Δploc⁡_r:
(1)Δpfull=Δpfric+Δploc⁡_r;Δpfric=λρ|w|22Dl;Δploc⁡_r=ζloc⁡_rρ|w|22,
where λ=λ(Re,Δ-) is the flow friction coefficient [2, 3]; Re=(ρ|w|D)/μ is the Reynolds criterion; ρ, w, and μ are the density, flow velocity, and dynamic viscosity coefficient of matter; D and l are the pipeline hydraulic diameter and length; Δ-=Δ/D is the relative pipe wall roughness; Δ is the absolute pipe wall roughness; and ζloc⁡_r is the local resistance coefficient. In a number of practically significant cases, the total resistance (total pressure losses) in calculations is treated as conventionally increased friction resistance [2]. In the first approximation, as applied to the problem of interest, let us use a similar approach; that is, we approximate the characteristic of interest (1) by the following relationship:
(2)Δpfull=λ(Re,δ-)ρ|w|22Dl,
where δ- is some increased (compared to Δ-) value of relative roughness that includes, among other things, the influence of local resistance. The parameter δ- is called equivalent relative roughness of the pipe wall. By analogy, the parameter δ=δ-D will be the equivalent absolute roughness of the pipe wall. Formula (2) has proven its credibility by many years of its application for the modeling of real gas mixture transport conditions in branched trunklines. This was primarily attributed to the actual closeness of the character of the above flow under design conditions to quadratic flow.Thus, the problem to be solved involves calculating the equivalent pipe wall roughnesses to incorporate the influence of flow resistance in segments of real gas networks [3, 5, 6]. The values of the above effective roughnesses are initially taken from the reference literature based on the available design and construction documentation. The values are then revised based on the consultations with the pipeline operating personnel and based on the results of a comparison of the calculated gas flow parameters (pressure, mass flow rate, and temperature) with corresponding field data gathered by remote measurement systems. It should be noted that the procedure of fitting the generalized semiempirical parameters of a gas pipeline system is basically performed at rather long time intervals—at most, once in two or three years—if there has been no repair or upgrade of the gas pipeline system over the given time period.In order to improve the quality and efficiency of the above fitting procedure and to minimize the influence of human intervention on the results of subsequent calculations, it makes sense to computerize it [4]. The basic prerequisite for such computerization as applied to a specific gas pipeline system is that the corresponding gas transport or gas distribution company should have an electronic data archive, including valve positions, gas withdrawal by consumers, and in-house gas consumption (here and below we do not differentiate between the terms “gas” and “gas mixture”), results of field measurements of basic gas flow parameters in the gas pipeline system, and parameters of gas transmission equipment. This data archive is generally collected and stored in the form of a computer database over a relatively long time period (one month to several years).


## 2. Preparation of Initial Data

To computerize the equivalent roughness fitting procedure, linear gas pipeline segments of the simulated gas pipeline network are conventionally divided inton nonoverlapping regions [4]. Neighboring of resulting pipeline regions is admissible (i.e., the regions can have a common boundary in space). The regions can be several hundred meters to several kilometers long. Then, for each resulting region of the linear gas pipeline segment, based on the reference literature, we establish the starting value of the equivalent roughness of the inner pipe wall.Based on the analysis of the electronic data archive, we identifyk nonoverlapping time intervals, over which gas transport was performed in company-specific non-steady-state nonisothermal regimes, which nevertheless differ from each other in their characteristics. We assume here that the gas pipeline system functioning over each of the identified time intervals was accompanied by a wide range of correct field measurements of gas flow parameters. The identification of the time intervals necessarily complies with the rule that the beginning of each identified time interval should be preceded by a period (though very small) corresponding to company-specific quasi-steady-state gas transport regimes. Violation of this rule will later make it impossible to use mathematical models of non-steady-state flow in gas pipeline systems to calculate target functions and functions in constraints of identification problems (see below). The above can be explained by the fact that, according to the laws of mathematical physics, solving non-steady-state problems requires setting not only boundary conditions but also initial conditions. In practice, initial data can be specified consistently only based on the results of preliminary solution of a corresponding steady-state problem.In our case, it is practical to identify time intervals such that they have the same lengthΔτ. If multiprocessor computers are used, this will enable efficient paralleling of the problem of gas dynamic model adjustment to real parameters of the gas pipeline system of interest. The recommended interval length Δτ can vary from several tens of minutes to several tens of hours. To improve the consistency of equivalent pipe wall roughness fitting, it makes sense to distribute the identified time intervals uniformly with respect to four calendar seasons. It is desirable for the time intervals to be spaced apart.As noted above, problem simulation by this method includes mathematical modeling of non-steady-state gas flow regimes for each of the identified time intervals with given values of equivalent pipe wall roughnesses in linear gas pipeline segments of the gas pipeline system of interest [1, 3]. Boundary conditions in modeling are specified at gas pipeline system boundaries in the form of various combinations of pressure, temperature, and gas flow rates.The modeling data are used to assess the closeness of the corresponding calculated and measured gas transport parameters in the gas pipeline system. Space-time pressure distributions we consider interest parameters. There are several reasons for this choice: the accuracy of pressure gauges is higher than the accuracy of flow rate and temperature meters [5, 6]; the number of pressure gauges in real networks is generally larger than the number of flow rate meters (as a rule, flow rate measurements are not taken inside the network); it follows from practical experience of extreme simulations that using a combination of pressure and flow rate parameters, which are fundamentally different in their scale, leads to significant worsening of the convergence of optimization simulations. The closeness of calculated and measured values is assessed for so-called identification points (IPs). IP locations are specified before equivalent roughness fitting based on the analysis of the layout of gas flow instruments in the pipeline network of interest and consistency and stability of their functioning. The preferred location of each IP should meet the following requirement [1, 3]: any significant change in time of rated gas dynamic regimes in the pipeline system of interest at the identification point must be accompanied by significant changes in time of gas flow parameters actually measured at this point. The distribution of IPs in the simulated pipeline network diagram should be as uniform as possible. Here we also assume that in the process of automatic adjustment of gas pipeline network models to real parameters of simulated objects the number of properly functioning IPs and their locations is fixed over all the identified time intervals.


## 3. Mathematical Formalization of Identification Problems

By analogy with [1], the closeness of calculated and measured time histories of the parameters under comparison at each IP will be analyzed in three senses [7]: closeness of the quality of two functional relationships (first-sense closeness); closeness of two functional relationships in the time-weighted average metric defined using the octahedral ∥⋯∥1 or Euclidean ∥⋯∥2 norm (second-sense closeness); closeness of two functional relationships within the framework of their uniform discrepancy, that is, in the metric ∥⋯∥0 (third-sense closeness).Calculated and measured time histories can be brought closer to each other in all the three senses by solving the classical identification problem:(3)∥pcalc(t,X)-pmeas⁡(t)∥ρ⟶min⁡X∈Ω⊂Rn,
where ∥⋯∥ρ is the vector norm, the type of which is determined by setting the parameter ρ=0,1,2; pcalc(t,X) is the vector function of calculated gas flow pressure values at the IP in the Euclidean space Rkm (it is practical to calculate these values using mathematical models presented in [1, 3]),
(4)pcalc(t,X)=([pcalc(t,X)]11,…,[pcalc(t,X)]k1,…,[pcalc(t,X)]km)T∈Rkm;pmeas⁡(t)=([pmeas⁡(t)]11,…,[pmeas⁡(t)]k1,…,[pmeas⁡(t)]km)T∈Rkm is the given vector function of measured gas flow pressure values at the IPs; t is the time (marching variable); X∈Ω⊂Rn is the vector of independent variables in the Euclidean space Rn (in the form of a set of sought equivalent inner wall roughnesses of pipeline segments); Ω={X∈Rn:a≤X≤b},  a∈Rn, b∈Rn are well-defined vectors establishing the boundaries in simple constraints on the range of admissible variations in the vector of independent variables; m is the number of defined IPs in the simulated pipeline system diagram; n is the number of independent variables.The problem (3) can have different statements depending on the chosen type of the vector function norm. For example, given the choice of the cubic norm ∥⋯∥0 we have a mini-max problem [8] with simple constrains on the variables:
(5)sup⁡(Δτ;1≤i≤k;1≤j≤m)|pcalc(t,X)-pmeas⁡(t)|ij⟶min⁡X∈Ω⊂Rn.
Solving the problem (5) provides the so-called uniform closeness of calculated gas transport parameters to their measured values. If we choose the octahedral vector norm ∥⋯∥1, the original problem (3) transforms into a nonlinear programming problem expressed as
(6)∫Δτ∑i=1k∑j=1m|pcalc(t,X)-pmeas⁡(t)|ijdt⟶min⁡X∈Ω⊂Rn.
Note that the problems (5) and (6) are equivalent only in terms of the definition of the global minimum. And the global solution of the problem (5) entails the global solution of (6). Unfortunately, we are not aware of scientifically based methods of search for the global extremum in such problems. An efficient search for the local solution to the problems (5) and (6) can be done by the known modified Lagrange functions (MLF) method supplemented with the variable metric method of search for the MLF minimum under simple constraints [9]. However, it is problematic to guarantee the coincidence of optimal (in the general case, local) solutions of the problems (5) and (6) in practice. The results of numerical solution of these problems can therefore differ from each other even if the same starting point is used.The choice of the Euclidean vector norm∥⋯∥2 in (3) results in the formulation of a new conditional optimization problem, which is practically equivalent to (6):
(7)∫Δτ∑i=1k∑j=1m[pcalc(t,X)-pmeas⁡(t)]ij2dt⟶min⁡X∈Ω⊂Rn.
Solving the problem (7) provides the root-mean-square closeness of calculated gas transport parameters to their measured values. It should be stressed here that the requirement of uniform closeness in (5) is stricter than the root-mean-square closeness in (7) [7]. As we show in [1, 3], it is also practical to solve the problem (7) by the MLF method.When solving the problem of fitting the equivalent roughnesses, one should give preference to the results of solving the problem (7) rather than the problem (5). This can be motivated by the fact that the compensation for individual spikes as a result of solving the problem (5) cannot guarantee satisfactory agreement between calculated and measured values in general. At the same time, partly smoothed individual discrepancy spikes in the functions under comparison that can theoretically remain after solving the problem (7) as a rule have a minor effect on the character of agreement between the functions of interest in general and can be ignored in the first approximation. Note that in most practical cases, given a fixed starting point, the runtime of the problem (7) (or (6)) is significantly smaller than that of the problem (5).Experience of production simulations involving the fitting of equivalent roughnesses of inner pipe walls shows that in the absence of sufficient IPs, in addition to the closeness of pressure values at the IPs, it makes sense to require that gas pressure differences along the pipelines between specifically identified pairs of IPs should also be close. In this case, the optimization problem will take the following form:(8)(∫Δτ∑i=1k∑r=1s[{pcalc(t,X)-pmeas⁡(t)}iMrpppppppppp+{pmeas⁡(t)-pcalc(t,X)}iNr]2dt(∫Δτ∑i=1k∑r=1s[{pcalc(t,X)-pmeas⁡(t)}iMr)1/2⟶min⁡X∈Ω⊂Rn,
where (Mr,Nr), r=1,s¯, is a predefined set of pairs of the numbers of IPs, between which the gas pressure difference along the pipeline is analyzed; s is the number of predefined pairs. The approach to this problem is similar to that of (7).Analysis of the closeness between calculated and measured time histories in the first sense (i.e., in essence, the closeness between the first derivatives of the functions under comparison) gives the following problem statement:(9)∫Δτ∑i=1k∑j=1m[∂pcalc(t,X)∂t-∂pmeas⁡(t)∂t]ij2dt⟶min⁡X∈Ω⊂Rn.


## 4. Analysis of Numerical Problem Solution

The values of equivalent inner wall roughness are determined by serial or parallel solution of the identification problems (7), (8), and (9) by the methods presented, for example, in [1, 3] as applied to pipeline transport. This gives three variants of calculated equivalent roughnesses. To choose the most acceptable set providing the highest possible accuracy of real gas flow identification by corresponding computational counterparts, one should use a criterion of comparison of different variants of computational gas dynamic flow regimes obtained for each of the k identified time intervals.Such a criterion was constructed based on the quantitative assessment in all the three senses of the closeness of calculated and measured time histories of gas pressure at each IP. The fundamentals and the order of its construction were presented in detail in [1]. It is therefore reasonable to restrict ourselves here only to the presentation of resulting expressions to calculate the generalized indicator of actual identification level by computational counterparts of real gas dynamic regimes of gas transport through the pipeline network of interest achieved in each identified time interval [1]:
(10)P_Identr=CGreenIPmin⁡_2mSgreenppppplppp×{Sgreen[P_Identgreen]r+Sblue[P_Identblue]rpppppplpppp+Sorange[P_Identorange]r},
where r=1,k¯. Consider(11)[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;[Ident_Level_2j]r=[Δτ-1∫Δτ|pcalc(t,X)-pmeas⁡(t)|jdt]r,j=1,m¯.[Ident_Level_2j]r, j=1,m¯, r=1,k¯, are the second-sense identification coefficients; Sgreen, Sblue, and Sorange are the scalar weight factors used to establish quantitative indicators of the identification level achieved, (Sgreen>Sblue≫Sorange>0); LGreenIPr is the number of IPs for the rth computational gas dynamic regime having the green identification level; KBlueIPr is the number of IPs for the rth computational gas dynamic regime having the blue identification level; NOrangeIPr is the number of IPs for the rth computational gas dynamic regime having the orange identification level. The conventional IP color code depending on the identification level achieved at the points was proposed in [1, 3]. According to this code, it is assumed that the identification level achieved in a small neighborhood of the jth IP for the time interval Δτ is high in the first and second senses if at least one of the following conditions is satisfied [1]: (1)
[Ident_Level_1j]r<CBlueIPmin⁡_1;(2)
[Ident_Level_2j]r≤CGreenIPmin⁡_2;(3)
[Ident_Level_1j]r≤CBlueIPmax⁡_1 and CGreenIPmin⁡_2 < [Ident_Level_2j]r < CBlueIPmin⁡_2,where CGreenIPmin⁡_2, CBlueIPmin⁡_1, CBlueIPmax⁡_1, and CBlueIPmin⁡_2 are specified empirical constants; [Ident_Level_1j]r, j=1,m¯, r=1,k¯, is the first-sense identification factor:
(12)[Ident_Level_1j]r=[Δτ-1∫Δτ|∂pcalc(t,X)∂t-∂pmeas⁡(t)∂t|jdt]r,pppppppppppppppppppj=1,m¯,r=1,k¯.
To render practical results in the IP diagram more informative, the point of interest is denoted by a green solid circle, and the identification level achieved at the point is called green. In a similar situation we assume that the identification level is satisfactory in the first and second senses if at least one of the following two composite conditions is satisfied (see (11) and (12)) [1]: (1)
CBlueIPmax⁡_1<[Ident_Level_1j]r and CGreenIPmin⁡_2<[Ident_Level_1j]r<CBlueIPmin⁡_2;(2)
CBlueIPmin⁡_1≤[Ident_Level_1j]r≤CBlueIPmax⁡_1 and CBlueIPmin⁡_2≤[Ident_Level_2j]r≤CBlueIPmax⁡_2,where CBlueIPmax⁡_2 is a specified empirical constant. In this case, the IP of interest in the IP layout diagram will be denoted by a blue solid circle (the blue identification level). Achievement of the identification level disputable in the first and second sense (orange identification level) is characterized by satisfying the following two conditions (see (11) and (12)) [1]: CBlueIPmin⁡_1≤[Ident_Level_1j]r≤CBlueIPmax⁡_1 and CBlueIPmax⁡_2<[Ident_Level_2j]r. If the above combinations of conditions are not fulfilled, we conclude that there is no identification in the first and second senses in the small neighborhood of the IP for the time interval Δτ. In the IP layout diagram, the IP is denoted by a red solid circle, and the absence of identification in its neighborhood corresponds to the red level.Analysis of the time history closeness in the first and second senses does not allow us to account for the influence of short-time sporadic spikes in measurement data that lead to discrepancies between compared time histories. Therefore, in the neighborhood of eachjth IP one should additionally analyze the closeness of the calculated and measured time histories in the third sense using a corresponding third-sense identification factor:
(13)[Ident_Level_3j]r={0,if[Ident_Level_1j]r≤CBlueIPmax⁡_1;sup⁡Δτ[|pcalc(t,X)-pmeas⁡(t)|j]r,otherwise,ppppppppppppppppppppj=1,m¯,r=1,k¯.
According to [1], the identification level established in the first and second senses should be downgraded: from green to blue if CGreenIPmin⁡_3<[Ident_Level_3j]r≤CIPsup⁡_3; from blue to orange if CBlueIPmin⁡_3<[Ident_Level_3j]r≤CIPsup⁡_3; from orange to red if COrangeIPmin⁡_3<[Ident_Level_3j]r≤CIPsup⁡_3; from any level to red if CIPsup⁡_3<[Ident_Level_3j]r, where CGreenIPmin⁡_3, CBlueIPmin⁡_3, COrangeIPmin⁡_3, and CIPsup⁡_3 are specified empirical constants. The downgrading of the identification level for each IP is a single step procedure; that is, the green level cannot be downgraded first to blue, then to orange, and finally to red.Thus, according to (10), for three variants of sets of equivalent inner wall roughnesses of pipeline segments, we calculate the values of the identification level indicator P_IdentrPROBLEM, where r=1,k¯; for PROBLEM=1, to calculate pressure values in (10), we used the set of roughness values obtained by solving (7); for PROBLEM = 2 we used values obtained by solving (8); and for PROBLEM = 3 we used values obtained by solving (9). It should be noted here that the identification level achieved grows with the value of the indicator P_IdentrPROBLEM. A comparative analysis of the roughness fitting results is done using arithmetic counters L, M, and N, corresponding to the sets of roughnesses obtained (initial values L:=0; M:=0; N:=0):
(14)L:=L+1,ifPIdentr1≥PIdentr2,P_Identr1≥P_Identr3;M:=M+1,ifPIdentr2≥PIdentr1,P_Identr2≥P_Identr3;N:=N+1,ifPIdentr3≥PIdentr1,lliP_Identr3≥P_Identr2,pppppppppppppppppppr=1,k¯.
As a solution to the problem of gas pipeline network model adaptation to real parameters of simulated objects in our case we choose the set of equivalent inner wall roughnesses of pipeline segments corresponding to the maximum value among the filled arithmetic counters L, M, and N. If the resulting parameters are L=0, M=0, and N=0, we conclude that the problem stated has no solution.


## 5. Results of Practical Application

To test the performance of the method proposed, we prepared sets of special test problems as applied to real gas pipeline layouts. Let us consider an example problem, namely, calculations of the parameters of the Moscow circular gas pipeline. The pipeline was conventionally divided into 26 segments of different lengths. The archive data bank was formed based on the pipeline operation data for three business days taken from different seasons of one year. The picture of pressures in the pipeline network at 98 IPS was produced by direct simulations for some roughness values, which were conventionally assumed to be “true.” For the model adaptation by the method proposed, initial values of roughness parameters were deliberately chosen to be different from “true” values. Following the calculations, we compared the calculated and the “true” values of equivalent roughness for each segment. The calculations demonstrated convergence of the method. The maximum relative error was 1.85%.The method of numerical automatic adaptation of constant characteristics of mathematical fluid transport models for trunkline and distribution pipeline systems has proven its efficiency as applied to production simulations for the Gazprom Regiongaz Moscow Limited Liability Company, where the method was used within the Alfargus knowledgeware package [10].


## 6. Conclusion

The paper describes a version of a practical method of mathematical gas flow model adaptation to actual parameters of real simulation objects. The method in itself is an extension of approaches to high-accuracy numerical recovery of transport flow in gas pipeline systems. This method can be used for the adaptation of constant characteristics of corresponding mathematical models to actual parameters of real engineering systems of the gas industry. This method has demonstrated its efficiency as applied to practical simulations by many years of its application to support pipeline operator decisions at gas transport and gas distribution companies and in the course of gas pipeline dispatcher training using research-grade computer simulators.



---

*Source: 146591-2014-02-09.xml*