Simulation and State Estimation of Transient Flow in Gas Pipeline

May 3, 2006 - Citation data is made available by participants in Crossref's Cited-by Linking service. For a more comprehensive list of citations to th...
0 downloads 16 Views 333KB Size
Download PDF
Ind. Eng. Chem. Res. 2006, 45, 3853-3863

3853

Simulation and State Estimation of Transient Flow in Gas Pipeline Networks Using a Transfer Function Model H. Prashanth Reddy,† Shankar Narasimhan,*,‡ and S. Murty Bhallamudi† Department of Chemical Engineering and Department of CiVil Engineering, Indian Institute of Technology, Madras, Chennai-600036, India

Dynamic simulation models of gas pipeline networks can be used for on-line applications such as state estimation, leak detection, etc. A prime requirement for such models is computational efficiency. In this paper, a transfer function model of a gas pipeline is used as a basis for developing a dynamic simulator for gas pipeline networks. The simulator is incorporated in a data reconciliation framework, which is ideally suited for on-line state estimation based on all available measurements of pressures and flow rates. The American Gas Association (AGA) model is used for making realistic computations of the gas compressibility. Accuracy and computational efficiency of the proposed method are evaluated by comparing our results with those obtained using a fully nonlinear second-order accurate finite difference method. The ability of the proposed approach for obtaining accurate state estimation from noisy measurements is demonstrated through simulations on an example network. We also demonstrate the use of the proposed approach for estimating an unknown demand at any node by exploiting the redundancy in measurements. 1. Introduction Pipeline networks are used extensively in all countries for transportation and distribution of natural gas and other light petroleum products for industrial and domestic use. These pipeline networks are large and intricate, consisting of pipes, compressors, valves, pressure reducers, flow limiters, etc. Transient flow conditions, which could typically last for several hours, are usually introduced in these systems during operation due to (i) changes in inlet and outlet flows/pressures, (ii) starting and stopping of compressors, and (iii) changes in control set points. Understanding such transient flow conditions through the use of dynamic models is therefore important for the design and safe operation of pipeline networks. Typically, in a dynamic simulation model, the governing continuity and momentum equations are solved by an appropriate numerical method. For a single pipe, these equations are a set of two nonlinear first-order hyperbolic partial differential equations, coupled with an equation of state for gas mixtures, and an equation for describing the friction loss. Equations of state could range from a simple ideal gas law to more elaborate models applicable to real gas mixtures.1 For a network, algebraic equations representing the compatibility conditions at the junctions and operating conditions of compressors and valves would also form part of the governing equations. Method of characteristics, explicit and implicit finite-difference methods, and finite-element methods may be used for solving these governing equations.2 Zhou and Adewumi3 presented hybrid total variation diminishing (TVD) schemes in which numerical oscillations resulting from dispersion errors are suppressed. Osiadacz4 compared different approximations of the governing equations and different numerical techniques for solving them. This study pertains to isothermal conditions, gases with a constant compressibility factor, and single pipelines. EmaraShabaik et al.5 also presented similar comparative investigations to evaluate the computational performance of CTCS (implicit central time and central space), CTAS (implicit central time * To whom correspondence should be addressed. E-mail: [email protected]. Phone: 91-44-22574165. † Department of Civil Engineering. ‡ Department of Chemical Engineering.

and alternating space), and MC (MacCormack predictor and corrector method) for their suitability for real time monitoring of fluid flow in pipelines, and they mentioned that the MC method gives more accurate results and faster dynamic response. Application of the method of characteristics becomes cumbersome for real gases because of the dependency of the compressibility factor on the pressure, temperature, and gas composition. Explicit finite difference methods require a small computational time step to ensure numerical stability. On the other hand, implicit finite difference methods involve inversion of large matrices during each time step, depending on the size of the network. These methods are very useful for accurate dynamic simulations and can be applied for off-line studies such as for design and safety analysis. However, they are not suitable for on-line monitoring applications such as state estimation and leak detection and identification because they are computationally very intensive. Kralik et al.6 developed a simplified transfer function model to investigate the errors introduced by the space discretization in finite-difference methods for solution of the governing differential equations. Results from the above study were later used by Kralik et al.7 to develop a computationally efficient implicit finite-difference method for simulation of large scale gas distribution systems. The sparsity of the matrix and structure of the network were exploited to further improve the computational efficiency. In our opinion, this is a significant work from the point of view of on-line implementations, although Kralik et al.7 themselves have used it only for off-line simulations. It may be noted that Kralik et al.7 have demonstrated the applicability of the method only for transients created by large changes in the boundary conditions such as the complete stoppage of inlet flow and rupture of the pipeline. Also, the compressibility factor was modeled as a simple linear function of the pressure. It may be noted that transfer function models have been developed in the past by other researchers also for analyzing fluid flow systems. Stecki and Davis8,9 reviewed and compared available linear partial differential equation (PDE) models for fluid transmission lines, and they have also derived a transfer function model relating fluid pressure P and volumetric flow rate Q at the two ends of a pipe. Makinen et al.10 presented a

10.1021/ie050755k CCC: $33.50 © 2006 American Chemical Society Published on Web 05/03/2006

3854

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006

similar rational transfer function model for dynamic simulation of pipelines transporting fluids. Matko et al.11 derived a nonlinear model in which a pipeline is represented by two transcendental transfer function equations relating flow variables at the inlet end to the outlet end. In all the above-mentioned models, density along the pipeline is taken to be constant and transfer function equations are used for the frequency domain analysis. In recent years, transfer function models have been developed and used for leak detection in water pipelines.12-17 All these methods employed a frequency domain analysis for identifying the leak. The objective of this study is to develop a method for state estimation in gas pipeline networks, based on available pressure and flow rate measurements, sampled at regular time intervals. For this purpose, a dynamic simulation model is combined with the concept of data reconciliation to exploit the redundancy in measured data.18 For lumped parameter systems which can be described by ordinary differential equations (ODEs), state estimation can be performed using the well-known Kalman filter (for linear systems) or extended Kalman filter (for nonlinear systems) techniques.19 To apply these methods for distributed systems such as gas pipelines, a spatial discretization procedure has to be used to convert the governing PDEs into ODEs. Such an approach has been used by Durgut and Leblebicioglu20 and Emara-Shabaik et al.21 for gas pipelines, and Estel et al.22 for heat exchangers. These methods typically result in high dimensional systems, which are not very suitable for on-line applications. To achieve computational efficiency, we propose using the approximate transfer function model developed by Kralik et al.6 for modeling the gas pipelines. The American Gas Association (AGA) model relating the compressibility factor to the pressure, temperature, and mixture composition is incorporated so that the proposed work has a wider applicability. The ultimate goal of our work is to use the proposed approach, which is based on linearization and approximate transfer functions, in an on-line model for state estimation during normal operating transient conditions and for detecting leaks caused by ruptures. These leaks do not introduce large or rapid changes in pressures and flows. The validity of the proposed approach for intended application is first assessed by comparing the results with those obtained using the second-order accurate McCormack explicit finite difference method.23 The ability of the proposed approach for obtaining accurate state estimation from noisy measurements is demonstrated through simulations on example networks. We also demonstrate the use of the proposed approach for estimating an unknown demand at any node by exploiting the redundancy in measurements. 2. Development of the Dynamic Simulation Model 2.1. Governing Equations for a Pipeline. The governing one-dimensional hyperbolic partial differential equations describing the conservation of mass and momentum for the transient flow of a gas through a constant diameter, rigid pipe are the following:2 Continuity Equation

∂p ∂M + (A/c2) ) 0 ∂x ∂t

(1)

Momentum Equation

1 ∂M ∂p + (g sin θ)p/c2 + λc2M|M|/(2DA2p) + )0 ∂x A ∂t

(2)

Equation of State

p/F ) zRT ) c2

(3)

where p is the pressure, M is the mass flow rate, g is the acceleration due to gravity, F is the density, θ is the inclination of the pipe, λ is the coefficient of friction, D is the inner diameter of the pipe, A is the cross-sectional area of the pipe, c is the pressure wave velocity, x is the distance along the pipeline, and t is time. The compressibility factor, z, is a function of pressure and temperature for a given gas. In this work, the compressibility factor is computed using the AGA1 model. The friction factor λ is determined using the Haaland explicit equation.24

6.9 /D 1.11 1 ) -1.8log10 + Re 3.7 xλ

[

( ) ]

(4)

where  is the roughness height of the pipe and Re is the Reynolds number. Equations 1-4 are standard equations, typically used for describing transients in natural gas pipeline networks.2 However, the following important points may be noted regarding the assumptions and limitations of the model. (i) Equation 4 is a steady-state friction model, which is adequate for our work, where the frequency in pressure and flow variations (due to normal operating transients and leaks) is not high. For rapid changes, frequency-dependent friction models as described by Zielke,25 Vardy and Brown,26 and Bergant et al.27 may be used. (ii) It may be noted that we are not using the energy equation for simulating the temperature variations. Instead, the effect of temperature variation along the pipe is taken into account through linear interpolation of available temperature measurements. It may be noted that this is not the same as assuming that the conditions are isothermal. Also, only within each sampling period, spatially averaged temperature for each pipe segment, based on currently measured values, is used while computing the compressibility. Typically, the temperature difference between the two end points of the pipeline is not very significant, especially for buried pipelines, in the case of normal operating transients. Therefore, the errors introduced due to this linearization procedure will not be significant. More discussion on this aspect can be found in Osiadacz and Chaczykowski.28 (iii) The flow velocity is assumed to be much smaller than the velocity of sound in the medium. Typically, flow velocities in gas pipelines are less than 10 m/s, while the velocity of sound is around 300 m/s. Equations 1-3 represent the description of the transient flow of a real gas in a single pipe. For a pipeline network, these equations along with (i) the specified boundary conditions, (ii) the continuity equation and pressure equality condition at junctions, and (iii) the equations describing the operating characteristics of nonpipe elements such as valves, compressors, etc. need to be solved simultaneously.2 Usually, this requires the conversion of eqs 1 and 2 to algebraic equations using appropriate spatial and temporal discretizations by a finitedifference approach. As noted earlier, finite-difference approaches require excessive computational time, which may prohibit the on-line deployment of such approaches. To overcome this difficulty, the transfer function model, developed by Kralik et al.,6 is used to represent the governing equations of a pipe. The transfer function model is based on linearization of the governing equations. Further more, first-order approximations of the transfer functions are required for obtaining time

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006 3855

domain solutions. Thus, this approach cannot provide as accurate results as the finite-difference based models. However, its main advantage is its computational efficiency, which makes it ideally suited for on-line applications. In fact, one of the objectives of this study is to get an idea about the accuracy of the transfer function based model by comparing the results with a benchmark finite-difference model. 2.2. Transfer Function Model. The transfer function model of a pipe relates the upstream pressure and mass flow rate to the downstream pressure and mass flow rate in the Laplace domain. Equations 1 and 2 are first expressed in terms of deviation variables from steady-state solution (∆M ) M - M0, ∆p ) p - p0). By assuming spatially averaged instantaneous values for the friction factor and the velocity of sound, the second and third terms in the momentum equation are linearized around the steady-state flow rate and average steady-state pressure in the pipe. This results in the following linear partial differential equations with time varying coefficients, similar to those obtained by Kralik et al.6

∂∆M ∂∆p + (AL/cj2) )0 ∂ξ ∂t

(5)

∂∆p + L(g sin θ)∆p/cj2 ∂ξ λLuj|uj| λLuj L ∂∆M ∆p + ∆M + ) 0 (6) 2 DA A ∂t 2Dcj h jc2 ) jzRT

(7)

uj ) M0cj02/(pj0A)

(8)

where ξ ) x/L denotes the nondimensional distance (L ) length of the pipe), the overbar denotes the spatially averaged values of the quantities in the pipe, and the subscript “0” denotes their values at the steady state. Taking the Laplace transform, we get the following two coupled linear ordinary differential equations.

[ ][

Equations 9 and 10 can now be written as

][

]

(11)

where

γ)

ALs λujL Ls -Lg sin θ λLuj|uj| + ; β) 2 + ; R) 2 2 DA A cj 2Dcj cj (12)

Equation 11 can be solved analytically to obtain the relations between the flows and pressures at the ends of the pipe in the Laplace domain.

∆M1 ) FM1,M2∆M2 + FM1,P1∆P1

FP2,P1 ) (e1/2γ2b)/(2b cosh b - γ sinh b); FP2,M2 ) (-2R sinh b)/(2b cosh b - γ sinh b) FM1,M2 ) (e -1/2γ2b)/(2b cosh b - γ sinh b); FM1,P1 ) (2β sinh b)/(2b cosh b - γ sinh b) (15) where

b ) (xγ2 + 4Rβ)/2

(13)

(16)

Although the above equations are derived to relate inlet and outlet conditions in a pipe, they can be used to relate conditions at any two positions along the pipe by substituting the appropriate value for the distance L. It is not possible to obtain analytical solutions in the time domain because it is difficult to obtain the inverse of the transfer functions in eq 15. Therefore, the transfer functions are approximated by a first-order series expansion of hyperbolic functions in eq 15 as follows.6

T1s 1 ; FM1,P1 ) ; FM1,M2 ) FP2,P1 ) k1 1 + Ts 1 + Ts (1 + T2s) 1 ; FP2,M2 ) -k2 (17) 1 + Ts 1 + Ts where

k 1 ) eΨ

(18)

λL 1 k2 ) eΨ/2 uj 1 + Ψ2 DA 24

(

)

(19)

1 λL2uj 1 1 - Ψ + Ψ2 T ) eΨ/2 2 6 24 2Dcj

(

(9)

d∆p(s,ξ) λLuj|uj| ∆p(s,ξ) + + L(g sin θ)∆p(s,ξ)/cj2 dξ 2Dcj2 λLuj L ∆M(s,ξ) + s∆M(s,ξ) ) 0 (10) DA A

(14)

where subscripts “1” and “2” denote the upstream and downstream ends of a pipe, respectively, and

T1 ) eΨ/2

d∆M(s,ξ) + (AL/cj2)s∆p(s,ξ) ) 0 dξ

d∆p(s,ξ) γ -R ∆p(s,ξ) ∂ξ ) d∆M(s,ξ) -β 0 ∆M(s,ξ) ∂ξ

∆P2 ) FP2,P1∆P1 + FP2,M2∆M2

T2 )

)

AL 1 1 + Ψ2 2 24 cj

(

(

1 λL2|uj| D + λ|uj| 6 Dcj2

Ψ)

)

(

1

(21)

))

Ψ2 1+ 24

λL uj|uj| gL sin θ 2D cj2 cj2

(20)

(22)

(23)

Kralik et al.6 have performed a frequency domain analysis to estimate the accuracy of the above simplified model in detail. They have recommended that the pipe segment length should be less than 45 000xD (length and diameter are in units of meters), which works out to be 25-45 km for typical diameters used in industry. This implies that long pipes should be divided into smaller segments for applying the above model. It should be noted that although Kralik et al.6 derived the above equations, they did not attempt to obtain the time domain solutions, which relate the upstream pressure and flow rate to the downstream pressure and flow rate for a pipe. Their objective was only to estimate the errors involved in the spatial discretization using frequency domain analysis. On the other hand, our objective is to use the above simplified model to solve the state estimation problem in the time domain. Therefore, eqs 17 are substituted in eqs 13 and 14, which are then analytically inverted. For this purpose, the time domain functions for ∆Ps and ∆Ms are first

3856

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006

2.3. Formulation of the State Estimation Problem. Equations 25 and 26 and the compatibility equations at the junctions for a network consisting of only pipes make up a linear system of equations. These can be expressed in matrix form as follows.

Am + Bu ) 0

Figure 1. Schematic of network.

expressed using their discrete values at regularly spaced time intervals, Ts. For example,

∆P1(t) ) ∆P1(kTs) kTs e t < (k + 1)Ts

(24)

The convolution theorem is then applied to obtain the Laplace inverse of eqs 14 and 13 to obtain the following equations in the time domain. N

∆p2(NTs) ) N

∑ i)1

[k1e-(N-i)T /T(1 - e-(T /T))∆p1(iTs)] + ∑ i)1 s

s

[

-k2T2 (1-e-(Ts/T)) T N k2T2 [ ∆M2(iTs)e-(N-i)Ts/T] + ∆M2(NTs) (25) T i)1

[-k2e-(N-i)Ts/T(1-e-(Ts/T))∆M2(iTs)]-

]

∑

(27)

In eq 27, the vector m consists of all measured variables (corresponding to time instants (N - n)Ts to NTs, where n ) number of past data samples required) and the vector u consists of all unmeasured variables (corresponding to time instants (N - n)Ts to NTs). The matrices A and B depend on the pipe parameters, sampling period, compressibility factor, and friction factor. Matrices A and B are time-dependent because the compressibility and friction factor vary with pressure, which is time-dependent. A simple network is considered in Appendix B to illustrate how eq 27 is obtained for a network using the transfer function model. In the framework of the state estimation problem, the adopted definition of the best estimate of the state is that which minimizes the quadratic difference between the measured values and the estimated values. In this context, the best estimates of variables m and u, for a given set of noisy measurements y for the variables m, can be obtained by solving the following weighted least squares estimation problem.

min (y - m)TQ-1(y - m) m,u

subject to

N

∆M1(NTs) )

[

[e-(N-i)T /T(1 - e-T /T)∆M2(iTs)] + ∑ i)1 s

]

N T1 T1 - (1 - e-Ts/T)[ ∆p1(iTs)e-(N-i)Ts/T] + ∆p1(NTs) T T i)1 (26)

∑

Am + Bu ) 0

s

where, current time t ) N × Ts. It is to be noted that the summation terms on the right-hand side of the above equations are weighted sums of the discrete values of pressures and flows at all times from the initial to the current time. However, the weights decrease exponentially from the current time, and thus, we need to consider only a finite number of time instances from the past (approximately 5 × T/Ts). For a typical pipe segment length of 10 km, the number of past samples required is about 200 when the sampling interval is 1 s. The complete discrete model in the time domain for the entire network is obtained by combining eqs 25 and 26 for all the pipe elements with (i) the continuity equation and pressure equilibrium equations at the junctions, (ii) the continuity equation and pressure drop equations at valves, (iii) the continuity equation and equations describing compressor operation, and (iv) the boundary conditions. As mentioned earlier, the governing partial differential equations are hyperbolic in nature. The characteristic theory23 shows that two boundary conditions should be specified for a well-posed problem for a single pipeline. For the case of subsonic flows, such as those considered in the present study, one boundary condition should be specified at the upstream end and one boundary condition should be specified at the downstream end. In the case of a network, one boundary condition should be specified at every inflow and outflow point. Also, demands at every internal node should be specified if there are withdrawals at these nodes. Thus, for the case of the network shown in Figure 1, three conditions, one at node 1, one at node 8, and one at node 5, should be specified.

where the matrix Q is the covariance matrix of errors in measurement. In the simplest case, this matrix is a diagonal matrix, with the diagonal element being the variance of the error of the corresponding measurement. The above problem is a standard linear estimation (reconciliation) problem for which the solution can be obtained by any nonlinear optimization technique. The dimensionality of the above optimization problem depends on the number of past time instances, n, and the size of the network. For on-line state estimation, the above optimization problem has to be set up and solved for each time instant in a moving horizon framework. This methodology would be computationally very intensive. Therefore, we obtain the state estimate for the current time instant NTs in a recursive manner by utilizing the estimates obtained for all the previous times. The optimization problem can be recast as given below.

j )TQ h -1(yj - m j) min (yj - m m j ,uj

A hm j +B h uj ) c

(28)

where m j and uj are the measured and unmeasured variables corresponding to the current time instant NTs, yj is the subvector of y corresponding to the measurements at the current time instant. A h, B h , and Q h are the corresponding submatrices of A, B, and Q, respectively. Vector c contains the weighted sum of the estimated flows and pressures at previous times. The simple network example is again used to illustrate how eq 28 is obtained from eq 27, as presented in Appendix B. An efficient approach to solve the above optimization problem is to make use of Crowe’s projection matrix technique as discussed by Narasimhan and Jordache.18 In this technique, the variables uj are first eliminated by constructing a projection h ) 0, where matrix P of dimension (ne - nu) × ne such that PB

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006 3857

ne is the number of equations and nu is the number of unmeasured variables in eq 28. The projection matrix P is obtained by singular value decomposition of the B h matrix. Thus, the above optimization problem becomes m j

subject to

(29)

The optimal estimates for the measured variables can be analytically obtained as

m jˆ ) yj - Q h (PA h )T[(PA h )Q h (PA h )T]-1P(A h yj - c)

available measurements at nodes run

pressure flow rate demand

h )Q h (PA h )T]-1PA h N)I-Q h (PA h )T[(PA

(32)

If just enough variables are specified, then the solution to the above problem corresponds to an objective function value of zero. In other words, we need to solve eqs 28 by substituting the specified values of m j . This would correspond to a simulation problem. If more measurements (specifications) than the minimum required to solve the problem are given, then the above formulation gives a best fit solution, taking into account the inaccuracies in the measured data. This would correspond to a state estimation problem. The boundary conditions for dynamic simulations described in section 2.2 translate into the observability conditions, i.e., minimum number and location of measurements for solving the state estimation problem. Thus, for a pipeline system such as that shown in Figure 1, a minimum of three measurements, one at node 1, one at node 5, and one at node 8, are required. However, in the state estimation framework, if the measurement for demand at node 5 is not available, then it is still possible to estimate both the state and the unknown demand, provided at least one additional measurement at any location in the network is available. In this regard, the state estimation problem is equivalent to the parameter estimation problem. The above discussion for observability conditions can be generalized as follows: (i) at least one measurement should be provided at each of the inlet and outlet boundaries; (ii) as many measurements as the number of internal demand points should be provided, which could be located at any point in the network. 3. Results and Discussion Accuracy and applicability of the proposed approach is evaluated through simulations on the hypothetical network shown in Figure 1. The hypothetical network (Figure 1) consists of 8 nodes including one source node, two delivery nodes, and 9 pipes. Tables 2 and 3 present the details of nodes and pipes, respectively. Natural gas (composition by mole fraction %: methane 93.42, nitrogen 0.12, carbon dioxide 2.36, ethane 1.76, propane 1.35, i-butane 0.31, n-butane 0.32, i-pentane 0.01, n-pentane 0.08, and n-hexane 0.01) having a viscosity of 0.000 01 N s/m2 is considered. The AGA gas model is used to

1 1

demand (MMSCMD)

elevation (m) MSL

pressure (kg/cm2)

temperature (K)

1 2 3 4 5 6 7 8

0 0 0 0 0.501 12 0 0 1.022 4

35 36 37 37 38 38 39 40

60.000 59.736 59.432 59.363 59.309 59.323 59.282 59.125

300.0 300.0 300.0 300.0 300.0 300.0 300.0 300.0

Table 3. Flow Rates throughout the Network in Steady-State Simulation pipe from to flow rate no. node no. node no. (MMSCMD)

where

5, 8 5, 8 5, 8 8

node no.

(30)

(31)

1 1 1, 3, 5, 7 1, 3, 5, 7

Table 2. Pressures at Nodes in Steady-State Simulation

The optimal estimates for the unmeasured variables uj can be obtained by substituting the estimates of the measured variables given by eq 30 in eq 28. The covariance matrix of errors in the estimates of measured variables can be obtained as follows:

h NT E[(m jˆ - m j )(m jˆ - m j )T] ) NQ

purpose

case 1 dynamic simulation case 2 state estimation nonredundant redundant case 3 unknown demand estimation

min (yj - m j )TQ h -1(yj - m j)

PA hm j ) Pc

Table 1. Details of Different Simulation Runs

1 2 3 4 5 6 7 8 9

1 2 3 4 5 3 5 6 7

2 3 4 5 6 6 7 7 8

1.5235 1.5235 0.6723 0.6723 -0.3052 0.8512 0.4764 0.5459 1.0224

length (m)

diam (m)

roughness height (m)

8 000.0 9 000.0 10 000.0 8 000.0 9 000.0 10 000.0 8 000.0 9 000.0 10 000.0

0.406 0.406 0.406 0.406 0.406 0.406 0.406 0.406 0.406

0.000 046 0.000 046 0.000 046 0.000 046 0.000 046 0.000 046 0.000 046 0.000 046 0.000 046

compute the compressibility factor. A constant ambient temperature of 300 K is specified at all points in the network. Details of different simulation runs, which are made on the above network, are presented in Table 1 for ready reference. 3.1. Case 1: Dynamic Simulation. We first test the accuracy of the proposed methodology by comparing the results for the case of a transient flow in the example network, with the results obtained using a second-order explicit finite-difference method. This finite-difference approach, known as the MacCormack method,23 is a classical technique for solving any nonlinear hyperbolic system of partial differential equations and is, therefore, chosen as a benchmark. The finite-difference solutions are obtained using a time step of 1 s and a spatial discretization of 375 m, which satisfy the requirement for numerical stability and grid convergence. In this simulation, the network (Figure 1) is assumed to be initially at a steady state corresponding to a specified upstream pressure of 60 kg/cm2 at node 1 and constant demands of 0.501 12 MMSCMD (million metric standard cubic meters per day) and 1.0224 MMSCMD at delivery nodes 5 and 8, respectively. Tables 2 and 3 present the initial steady-state pressures and flows obtained for the network. The negative flow in the pipe element number 5 implies that the flow direction is from node 6 to node 5. In the present case, starting from the above-mentioned steady state, the demand at node 8 varies as shown in Figure 2, while the pressure at node 1 and the demand at node 5 are constant. These three variables are treated as measured (without any noise) and are provided as inputs to the estimator. Since these constitute the minimum required for state estimation, this case corresponds to the standard dynamic simulation of the network. The dynamic simulations are carried out for 6000 s.

3858

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006

Figure 2. Specified demand variation at node 8. Figure 5. Comparison of pressure at node 5, case 1.

Figure 3. Comparison of mass flow rate at node 1, case 1. Figure 6. Comparison of pressure at node 8, case 1.

Figure 4. Comparison of pressure at node 2, case 1.

Simulations with the proposed method are carried out for two sampling intervals of 1 s and 15 s in order to demonstrate the effect of sampling frequency on the accuracy of the results. Estimated pressures at nodes 2, 5, and 8 and the estimated mass flow rate at node 1 obtained using the proposed method are compared with the benchmark solution in Figures 3-6, respectively. It is observed from these figures that the proposed approach simulates the dynamic conditions in the network accurately, when the sampling interval is 1 s. However, there is an error of approximately 4% in the peak mass flow rate value when the sampling interval is increased to 15 s. The error in the simulated peak pressure values is approximately 0.1%. It is also noted that there is a phase lag between the results obtained by the proposed model and the benchmark solution, when the sampling interval is 15 s. These phase and amplitude errors are due to the linearization of the governing equations and first-order approximations of the transfer functions. Despite the errors, the trends are well simulated by the proposed approach. Numerical experimentation showed that the results

obtained using a sampling interval of 5 s are comparable to those obtained using a sampling interval of 1 s. It should be noted that, in the proposed approach, the pipe elements are not discretized, while, in the approach using the finite-difference technique, the pipes are discretized using a space step of 375 m. A comparison of the computational time requirement on an Intel Pentium IV 3.06 Hz processor shows that the proposed approach consumes 34 min, whereas the benchmark solution consumes 14 h, for a sampling interval of 1 s and a total simulation time of 6000 s. In other words, the proposed approach requires only 0.34 s to process data with a 1-s sampling interval, whereas the benchmark solution has taken 8.4 s. Thus, the proposed approach can be used in on-line applications, for small sampling intervals without compromising accuracy. A sensitivity analysis (results not shown here) indicates that there is no need to subdivide even pipes of 30 km in length. However, as already pointed out, the sampling frequency should be sufficiently high in order to obtain accurate results. Studies conducted using the ideal gas model instead of the AGA model also led to the same conclusions. 3.2. Case 2: State Estimation. In the preceding subsection, we demonstrated the accuracy of the proposed methodology for simulating transient conditions in a network. One of the main advantages of our approach is its ability to consider all available measured data in order to obtain an optimal state estimate. Typically, measurements contain errors, and our approach uses the redundancy in data availability to reduce the effect of such errors in state estimation. This case study illustrates the applicability of our approach for the above-mentioned purpose. The example network shown in Figure 1 is considered. To generate synthetic measured data, transient simulations are carried out, as described in the previous section, for the demand transient at node 8 (Figure 2). This gives the “true” flows and

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006 3859

Figure 7. True and estimated pressures at node 1, case 2.

Figure 8. True and estimated mass flow rates at node 1, case 2.

Table 4. Reduction in RMS Error with Increased Redundancy variable

node

just specified with noise

redundant measurements with noise

pressure flow pressure flow pressure pressure flow pressure demand pressure pressure pressure demand

1 1 2 2 3 4 4 5 5 6 7 8 8

0.010 0.055 0.025 0.005 0.025 0.024 0.001 0.024 0.003 0.024 0.024 0.021 0.003

0.005 0.000 0.019 0.000 0.018 0.017 0.0001 0.017 0.001 0.018 0.017 0.015 0.001

pressures at all nodes, at a sampling interval of 1 s. The “measurements” corresponding to any specified set of measured variables are obtained by adding “noise” to the above-simulated true values. The noises are generated randomly from the Gaussian distribution with zero means and standard deviations of 10-2 kg/cm2 and 3.0 × 10-3 MMSCMD in pressure and flow rate, respectively. In the nonredundant case, noisy measurements of pressure at node 1 and demands at nodes 5 and 8 are provided as inputs to our approach and the flows and pressures are estimated in the network. On the other hand, in the redundant case, noisy measurements of pressures at nodes 1, 3, 5, and 7, inflow at node 1, and demands at nodes 5 and 8 are provided as input. Estimated flows and pressures from the above two runs are compared with the true values in order to evaluate the data reconciliation ability of the proposed approach. Significant improvement in the estimates of all flows and pressures is obtained in the redundant case as indicated by the root mean squared (RMS) error values presented in Table 4. Figures 7, 8, 9, 10, and 11 show the comparison between the true values and estimated values for the redundant and nonredundant cases, for pressure at node 1, flow rate at node 1, flow rate at node 2, demand at node 5, and demand at node 8, respectively. It is very clear that the proposed methodology is able to use all the available measurements to accurately estimate the state of the system, even when the measurements are noisy. However, a closer inspection of the results reveals a small, but negligible, bias between the estimated and true quantities. This may be attributed to (i) approximations involved in linearizing eqs 1 and 2 and (ii) first-order approximations of the transfer functions in eq 15. 3.3. Case 3: Estimation of Unknown Demand. Typically, in any pipeline system, the flow rates and pressures at all supply

Figure 9. Estimated mass flow rates at node 2 with and without redundancy in measurements, case 2.

Figure 10. True and estimated demands at node 5, case 2.

and delivery points are measured. However, in case there is a leak, the outflow due to the leak is unknown. The state estimation approach can be used to estimate the leak magnitude by treating it as an unknown demand, as long as additional measurements of pressures are available along the pipe. This capability of the proposed state estimation approach is illustrated below. For simplicity, the location where the leak (unknown demand) occurs is assumed to be known. Noisy measurements as in case 2 are generated with the exception that a step change in demand at node 5 from 0.5011 to 0.5512 MMSCMD is introduced at time t ) 1000 s. The above-mentioned demand is assumed to be unmeasured and is estimated from the other measurements using the proposed

3860

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006

accurately if there are no transients induced by demand variations at other nodes. The above-mentioned capability of the proposed method is important for extending the approach to leak detection and isolation in gas pipeline networks, wherein even the location of the leak has to be estimated along with the magnitude. We have recently developed such an approach, the results of which will be communicated in future. 4. Conclusions

Figure 11. Estimated demands at node 8 with and without redundancy in measurements, case 2.

A transfer function model for a single gas pipeline has been successfully used as a basis for developing an accurate and computationally efficient approach for dynamic simulations in gas pipeline networks. The simulator has been incorporated into a data reconciliation framework for online state estimation. Test problems on an example network indicated that the proposed method was 25 times faster than the explicit finite-difference approach. We also demonstrated that the proposed approach can be used to estimate unknown demands. The abovementioned features, coupled with the computational efficiency, make the approach ideally suited for on-line leak detection and identification. Acknowledgment This research work was financially supported by GAIL (India) Ltd. under the sponsored project “Development of leak detection methods in gas pipeline networks”. Appendix A: Derivation of Equations 25 and 26

Figure 12. Estimated demand at node 5, case 3 with demand transient at node 8.

Steps followed to invert eqs 13 and 14 to obtain eqs 25 and 26 are presented here. Consider eq 14 as shown below.

∆P2 ) FP2,P1∆P1 + FP2,M2∆M2

(A-1)

where

(1 + T2s) 1 FP2,P1 ) k1 ; FP2,M2 ) -k2 1 + Ts 1 + Ts

Figure 13. Estimated demand at node 5, case 3 with no demand transient.

approach as a part of state estimation described in the earlier section. The simulation is carried out for 6000 s, and the estimated demand at node 5 is presented in Figure 12. Prior to the step change, the estimated demand at node 5 is approximately 0.50 MMSCMD, which is close to the true value. Subsequently, the mean estimated demand differs from the true value due to the transient introduced in demand at node 8. However, the mean estimated demand eventually converges to the correct mean value of 0.55 MMSCMD. In another case, noisy measurements are generated for a step change in demand at node 5 as before, but for a constant demand of 1.022 MMSCMD at node 8; i.e., there is no demand variation at node 8. For this case, the unknown demand estimated at node 5 is presented in Figure 13. It can be seen from this figure that the mean estimated unknown demand is close to the true value at all times. It is clear from these two tests that the proposed approach is able to reconcile the unknown demand. It is also noted that the unknown demand at node 5 is estimated more

and s is the Laplace operator. In eq A-1, all the variables ∆P2, ∆P1, ∆M2, FP2,P1, and FP2,M2 are in the Laplace domain. Therefore, a convolution integral is required to obtain the Laplace inverse of both the terms on the RHS of eq A-1. For the sake of brevity, only the Laplace inverse of the first term is presented here. A similar procedure is adopted to obtain the Laplace inverse of other terms in this as well as in eq 13. The convolution theorem states that

L-1{F(s) G(s)} )

∫0tf(t - τ) g(τ) dτ

(A-2)

Application of (A-2) to L-1{FP2,P1(s)∆P1(s)} and writing it for the current time instant t ) NTs gives

L-1{FP2,P1(s)∆P1(s)} )

k

∫0NT T1e-(NT -τ)/T∆P1(τ) dτ s

s

(A-3)

where

L-1

{

}

k1 k1 ) e-t/T 1 + Ts T

(A-4)

Note that ∆P and ∆M variables for time t ) 0 are zero because

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006 3861

following: 6 equations for pipe elements, 4 equations for two intermediate junction nodes, and 1 equation at the end node (the total number of equations at any time step ) 11). In eqs 25 and 26, t ) NTs represents the current time instant. Let us consider that N ) 3, for the sake of illustration. Equations at the current time instant, i.e., N ) 3, can be written as

Figure B1. Definition sketch for series pipeline.

Figure B2. Discretization of pipeline system shown in Figure B1.

we start from a steady state. The time domain functions for ∆P and ∆M are then approximated using their discrete values at regularly spaced time intervals, Ts. For example,

∆P1(t) ) ∆P1(kTs) kTs e t < (k + 1)Ts The RHS of eq A-3 becomes

{

∫

∫

Ts

2Ts

-(NTs-τ)/T ∆P1(Ts) dτ + T e-(NTs-τ)/T∆P1 k1 0 e s ) T (2T ) dτ + ... + NTs e-(NTs-τ)/T∆P (NT ) dτ s 1 s (N-1)T

k1 ) T

∫

{

s

(A-5)

}

T(e-((N-1)Ts/T) - e-(NTs/T))∆P1Ts + T(e -((N-2)Ts/T) e-((N-1)Ts/T))∆P12Ts + ... ... + T(1 - e-(Ts/T))∆P1(NTs)

) k1{e-(N-1)Ts/T(1 - e-(Ts/T))∆P1(Ts) + e-(N-2)Ts/T

L-1

{

k1

I101∆M2(2Ts) + I111∆P1(2Ts) + I102∆M2(Ts) + I112∆P1(Ts) ∆P2(3Ts) ) I120∆M2(3Ts) + I130∆P1(3Ts) + I121∆M2(2Ts) + I131∆P1(2Ts) + I122∆M2(Ts) + I132∆P1(Ts) ∆M3(3Ts) ) I200∆M4(3Ts) + I210∆P3(3Ts) + I201∆M4(2Ts) + I211∆P3(2Ts) + I202∆M4(Ts) + I212∆P3(Ts)

}

(1 - e-(Ts/T))∆P1(2Ts) + ... + T(1 - e-(Ts/T))∆P1(NTs)} Therefore,

∆M1(3Ts) ) I100∆M2(3Ts) + I110∆P1(3Ts) +

∆P4(3Ts) ) I220∆M4(3Ts) + I230∆P3(3Ts) + I221∆M4(2Ts) + I231∆P3(2Ts) + I222∆M4(Ts) + I232∆P3(Ts) ∆M5(3Ts) ) I300∆M6(3Ts) + I310∆P5(3Ts) + I301∆M6(2Ts) + I311∆P5(2Ts) + I302∆M6(Ts) + I312∆P5(Ts) ∆P6(3Ts) ) I320∆M6(3Ts) + I330∆P5(3Ts) + I321∆M6(2Ts) + I331∆P5(2Ts) + I322∆M6(Ts) + I332∆P5(Ts)

}

∆M2 - ∆M3 ) 0

∆P1(s) )

1 + Ts

i)N

k1{e ∑ i)1

-(((N-i)Ts)/T)

∆P2 - ∆P3 ) 0 (1 - e

-(Ts/T)

)∆P1(iTs)} (A-6)

∆M4 - ∆M5 ) 0 ∆P4 - ∆P5 ) 0

Similarly, the time domain functions can be derived for other terms in eqs 13 and 14. Appendix B: Derivation of Equations 27 and 28 Derivation of Equation 27. Explanation is given here on how to set up eqs 25 and 26 and compatibility conditions at the junctions together in the form of eq 27. A simple series pipeline consisting of three pipe elements (Figure B1) is considered for illustration. 1, 2, 3, and 4 are the node numbers. The pipe numbers are represented in parentheses. In this series pipeline, we have assumed that pressure measurements at nodes 1, 2, and 4, the mass flow rate measurement at node 1, and the demand measurement at node 4 are available. In the implementation of the methodology, the abovementioned example network is divided into individual pipe elements as shown in Figure B2. In this representation, nodes from 1 to 5 are associated with two variables (pressure and mass flow rate) and node 6 is associated with pressure, mass flow rate, and demand. Thus, the variables are the following: P1, M1, P2, M2, P3, M3, P4, M4, P5, M5, P6, M6, and D6 (P indicates pressure, M indicates mass flow rate, D indicates demand, and the subscript indicates the node number). Among these variables, P1, M1, P2, P6, and D6 are measured variables. At any time instant, two transfer function equations (in the time domain) are available for each pipe element, two equations (mass flow rate conservation and pressure equality) are available at each intermediate junction node, and one equation (mass flow rate conservation) is available at the end node. Therefore, equations available at any time instant for the system shown in Figure B2 are the

∆M6 - ∆D6 ) 0

(B-1)

Where

(

Ip0j ) 1.0 - exp

( )((

Ip1j ) -

T1 T

( )) ( ) ( )) ( )) Ts T

1.0 - exp

Ts exp -j T

Ts T

Ts exp -j T

and if j ) 0

Ip1j ) -

( )(( ( )) ( )) ( ) ( )( ( )) ( ) T1 T

1.0 - exp

Ip2j ) -k2 1.0 -

Ts T

Ts exp -j T

T2 Ts 1.0 - exp T T

+

exp -j

T1 T

Ts T

if j ) 0

(

Ip2j ) -k2 1.0 -

)(

( )) ( ) ( )) ( )

T2 Ts 1.0 - exp T T

(

Ip3j ) k1 1.0 - exp

Ts T

Ts T2 exp -j - k2 T T

Ts exp -j T

In the coefficient Ipqj, p denotes the corresponding pipe element number, q ) 0 and 1 denote the coefficient for mass flow rate and the coefficient for pressure, respectively, in eq 25, and q )

3862

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006

[ ] [ ] [ ] [] [] [ ] [ ]

2 and 3 denote the coefficient for mass flow rate and the coefficient for pressure, respectively, in eq 26. The superscript j denotes (N - i) in eqs 25 and 26 in which i is the coefficient at the ith time instant. Equations B-1 can be written for all the time instants 1, 2, and 3 and can be grouped together to obtain a matrix equation in the form of eq 27. For N ) 3 and the series pipeline shown in Figure B2, the matrix size of A is the following: number of rows ) 11 × 3 ) 33 (no. of equations at any time instant ) 11 and no. of time instances ) 3) and number of columns ) 5 × 3 ) 15 (no. of measured variables ) 5). The size of the B matrix is the following: number of rows ) 11 × 3 ) 33 and number of columns ) 8 × 3 ) 24 (no. of unmeasured variables ) 8). The size of vector m is 5 × 3 ) 15, and it contains the variables ∆P1, ∆M1, ∆P2, ∆P6, and ∆D6 at the three time instances 3Ts, 2Ts, and Ts. Correspondingly, the size of the vector u ) 8 × 3 ) 24. Derivation of Equation 28. It is assumed that the values of ∆P and ∆M at time instances 2Ts and Ts are known while writing eq B-1 for time instance 3Ts. These values are known from previous estimates. Equation B-1 at the current time can now be written in matrix form as shown below. 1.0-I110-I100 0 0 0 0 0 0 0 0 0 0 0 -I130-I1201.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 -I21 -I20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -I230-I2201.0 0 0 0 0 0 0 0 0 0 0 1.0 -I31 -I300 0 0 0 0 0 0 0 0 0 0 0 -I330-I3201.0 0 0 0 1.0 0 -1.0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 -1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 -1.0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 -1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 -1.0 ∆M1 ∆P1 ∆M2 ∆P2 ∆M3 ∆P3 ∆M4 ) ∆P4 ∆M5 ∆P5 ∆M6 ∆P6 ∆D6

I101∆M2(2Ts) + I111∆P1(2Ts) + I102∆M2(Ts) + I112∆P1(Ts) I121∆M2(2Ts) + I131∆P1(2Ts) + I122∆M2(Ts) + I132∆P1(Ts) I201∆M4(2Ts) + I211∆P3(2Ts) + I202∆M4(Ts) + I212∆P3(Ts) I221∆M4(2Ts) + I231∆P3(2Ts) + I222∆M4(Ts) + I232∆P3(Ts) I301∆M6(2Ts) + I311∆P5(2Ts) + I302∆M6(Ts) + I312∆P5(Ts) I321∆M6(2Ts) + I331∆P5(2Ts) + I322∆M6(Ts) + I332∆P5(Ts) 0 0 0 0 0

(B-2)

Equation B-2 can be written in the form of eq 28 by separating the corresponding matrices for measured and unmeasured variables. The resulting equation in the form of eq 28 can be written as shown below:

1.0 -I110 0 -I130 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 1.0 0 0 0 0 0 1.0 0 0 0

0 0 0 0 0 0 ∆M1 0 0 ∆P1 0 0 ∆P2 + 1.0 0 ∆P6 0 0 ∆D 6 0 0 0 0 0 0 0 -1.0

[]

-I100 0 0 0 0 -I12 0 0 0 0 1.0 -I210 -I200 0 0 -I230 -I220 0 0 0 0 0 0 0 0 1.0 -1.0 0 0 0 0 -1.0 0 0 0 0 1.0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 1.0 -I31 -I300 0 0 -I330 -I320 0 0 0 0 0 0 0 0 0 -1.0 0 0 1.0 0 -1.0 0 0 0 0 1.0

∆M2 ∆M3 ∆P3 ∆M4 ∆P4 ) ∆M5 ∆P5 ∆M6

I101∆M2(2Ts) + I111∆P1(2Ts) + I102∆M2(Ts) + I112∆P1(Ts) I121∆M2(2Ts) + I131∆P1(2Ts) + I122∆M2(Ts) + I132∆P1(Ts) I201∆M4(2Ts) + I211∆P3(2Ts) + I202∆M4(Ts) + I212∆P3(Ts) I221∆M4(2Ts) + I231∆P3(2Ts) + I222∆M4(Ts) + I232∆P3(Ts) I301∆M6(2Ts) + I311∆P5(2Ts) + I302∆M6(Ts) + I312∆P5(Ts) I321∆M6(2Ts) + I331∆P5(2Ts) + I322∆M6(Ts) + I332∆P5(Ts) 0 0 0 0 0

(B-3)

In the above equation, ∆P1, ∆M1, ∆P2, ∆P6, and ∆D6 are the measured variables and so appear in the vector “m j ” vector. ∆M2, ∆M3, ∆P3, ∆M4, ∆P4, ∆M5, ∆P5, and ∆M6 are the unmeasured variables and so appear in the vector uj. The right-hand vector, “c” contains the weighted sum of ∆Ps and ∆Ms at the previous time instances 2Ts and Ts.

Ind. Eng. Chem. Res., Vol. 45, No. 11, 2006 3863

Literature Cited (1) Compressibility Factor of Natural Gas and Related Hydrocarbon Gases, AGA Report No. 8, American Gas Association: Arlington, VA, 1994. (2) Osiadacz, A. J. Simulation and Analysis of Gas Networks; Gulf Publishing Company: Houston, TX, 1989. (3) Zhou, J.; Adewumi, M. A. Simulation of transients in natural gas pipe lines using hybrid TVD shemes. Int. J. Numer. Methods Fluids 2000, 32 (4), 407-437. (4) Osiadacz, A. J. Different transient modelsslimitations, advantages and disadvantages. In Proceedings of the PSIGsThe 28th Annual Meeting, San Francisco, 1996. (5) Emara-Shabaik, H. E.; Khulief, Y. A.; Hussaini, I. Simulation of transient flow in pipelines for computer-based operations monitoring. Int. J. Numer. Methods Fluids 2004, 44 (3), 257-275. (6) Kralik, J.; Stiegler, P.; Vostry, Z.; Zavorka, J. Modeling the dynamics of flow in gas pipeline. IEEE Trans. Syst., Manage., Cybern. 1984, SMC14 (4), 586-596. (7) Kralik, J.; Stiegler, P.; Vostry, Z.; Zavorka, J. A Universal Dynamic Simulation Model of Gas Pipeline networks. IEEE Trans. Syst., Manage., Cybern. 1984, SMC-14 (4), 597-606. (8) Stecki, J.; Davis, D. Fluid Transmission Lines-Distributed Parameter Models, Part 1: A Review of the State of the Art. Proc. Inst. Mech. Eng., Part A 1986, 200, 215-228. (9) Stecki, J.; Davis, D. Fluid Transmission Lines-Distributed Parameter Models, Part 2: Comparison of Models. Proc. Inst. Mech. Eng., Part A 1986, 200, 229-236. (10) Makinen, J.; Piche, R.; Ellman, A. Fluid Transmission Line modelling Using a Variational method. Mathematics report 72; Tampere University of Technology: Tampere, Finland, 1996. (11) Matko, D.; Geiger, G.; Gregoritza, W. Pipeline simulation techniques. Math. Comput. Simul. 2000, 52 (3-4), 211-230. (12) Ferrante, M.; Brunone, B. Pipe system diagnosis and leak detection by unsteady-state tests. 1. Harmonic analysis. AdV. Water Resour. 2003, 26, 95-105. (13) Ferrante, M.; Brunone, B. Pipe system diagnosis and leak detection by unsteady-state tests. 2. Wavelet analysis. AdV. Water Resour. 2003, 26, 107-116. (14) Kim, S. H. Extensive Development of leak detection algorithm by impulse response method. J. Hydraulic Eng. 2005, 131 (3), 201208.

(15) Lee, P. J.; Vitkovsky, J. P.; Lambert, M. F.; Simpson, A. R.; Liggett, J. A. Frequency domain analysis for detecting pipeline leaks. J. Hydraulic Eng. 2005, 131 (7), 596-604. (16) Mpesha, W.; Gassman, S. L.; Chaudhry, M. H. Leak detection in pipes by frequency response method. J. Hydraulic Eng. 2001, 127 (2), 134147. (17) Mpesha, W.; Chaudhry, M. H.; Gassman, S. L. Leak detection in pipes by frequency response method using a step excitation. J. Hydraulic Res. 2002, 40 (1), 55-62. (18) Narasimhan, S.; Jordache. C. Data Reconciliation & Gross Error Detection - An Intelligent Use of Process Data; Gulf Publishing Company: Houston, TX, 2000. (19) Muske, K. R.; Edgar, T. F. Nonlinear state estimation. In Nonlinear Process Control; Henson, M. A., Seborg, D. E., Eds.; Prentice-Hall: New Jersey, 1997; pp 311-370. (20) Durgut, I.; Leblebicioglu, K. Kalman-filter-based observer design around optimal control policy for gas pipelines. Int. J. Numer. Methods Fluids 2000, 32 (4), 407-437. (21) Emara-Shabaik, H. E.; Khulief, Y. A.; Hussaini, I. A nonlinear multiple-model state estimation scheme for pipeline leak detection and isolation. Proc. Inst. Mech. Eng., Part I 2002, 216 (10), 497-512. (22) Estel, L.; Bagui, F.; Abdelghani-Idrissi, M. A.; Thenard, C. Distributed state estimation of a counter current heat exchanger under varying flow rate. Comput. Chem. Eng. 2000, 24 (1), 53-60. (23) Chung, T. J. Computational Fluid Dynamics; Cambridge University Press: Cambridge, 2002. (24) White, F. M. Fluid Mechanics, 2nd ed.; MacGraw-Hill Book Company: New York, 1986. (25) Zielke, W. Frequency-dependent friction in transient pipe flow. ASME J. Basic Eng. 1968, 90, 109-115. (26) Vardy, A. E.; Brown, J. M. B. Transient turbulent friction in fully rough pipe flows. J. Sound Vib. 2004, 270 (1-2), 233-257. (27) Bergant, A.; Simpson, A. R.; Vutkovsky, J. Developments in unsteady pipe flow friction modeling. J. Hydraulic Res. 2001, 39 (3), 249257. (28) Osiadacz, A. J.; Chaczykowski, M. Comparison of isothermal and nonisothermal pipeline gas flow models. Chem. Eng. J. 2001, 81, 4151.

ReceiVed for reView June 24, 2005 ReVised manuscript receiVed February 15, 2006 Accepted March 10, 2006 IE050755K