- Email: [email protected]
A linear-rate analog approach for the analysis of natural gas transportation networks
Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
A linear-rate analog approach for the analysis of natural gas transportation networks Qian Sun, Luis F. Ayala* John and Willie Leone Family Department of Energy and Mineral Engineering, The Pennsylvania State University, University Park, PA 16802, United States
a r t i c l e i n f o
a b s t r a c t
Article history: Received 26 October 2016 Received in revised form 18 January 2017 Accepted 27 March 2017 Available online 6 April 2017
Modeling natural gas transportation networks poses a number of challenges due to the significant nonlinearities associated to the governing equations. Pressure (nodal or p-) formulations and flow (nodalloop or q-) formulations are the most commonly deployed approaches used to formulate gas flow network models. They treat nodal pressures and pipe branch flow rates as primary unknowns, respectively, and the Newton-Raphson method is the typical choice used to solve the resulting system of pipe network equations. The major disadvantage of Newton-Raphson methods is their tendency to hopelessly diverge when good initializations for the unknown pressure and flow variables are not available. In order to overcome this drawback, a linear-pressure analog approach, applicable to the p-formulation, was recently proposed to formulate an initial-guess-free solution protocol. However, solving the q-formulation rather than the p-formulation would be an alternate, practical, and desirable way to study gas flow in pipeline networks because the system of equations is predominantly lineardwith the exception of the loop equations. Linear Theory and Hardy Cross methods have been used in the past solve such qformulation. Unfortunately, these methods are not initial-guess-free and have been shown to become potentially unstable and thus inefficient because their convergence is also strongly associated with the availability of good initial guesses. This work proposes a linear-rate analog method capable of effectively solving the q-formulation. Through case studies, the proposed linear-rate analog method is shown to be a robust and initial-guess-free solution scheme to effectively model horizontal and inclined natural gas pipeline networks. © 2017 Elsevier B.V. All rights reserved.
Keywords: q-formulation Gas flow in conduits Networks Linear-rate analog
1. Introduction Reliable modeling of natural gas transportation networks continues to increase in importance as a result of increased energy demands and production from large unconventional gas fields (EIA, 2013; Makogon et al., 2007; Wood et al., 2008). Solving a natural gas pipeline network system is challenging due to the high nonlinearity of the governing equations of gas flow in conduits. Therefore, the development of robust mathematical methods to model gas flow in pipeline becomes of utmost importance for the study of optimized design and performance. Governing equations of gas pipe flow in networks are derived from energy conservation statements for steady state conditions (Ayala, 2013). Typically, pipe networks consist of nodes and pipe branches. External demands and supplies can be present at the
* Corresponding author. E-mail address: [email protected] (L.F. Ayala). http://dx.doi.org/10.1016/j.jngse.2017.03.027 1875-5100/© 2017 Elsevier B.V. All rights reserved.
nodes. In these networks, closed circuits formed by pipe branches are known as pipe loops. The formulation of gas flow network equations can be categorized as pressure (p-), pipe flow (q-) or loop-rate correction (Dq-) approaches, where node pressures, branch flow rates or corrections to flow rates within a loop, respectively, are treated as primary unknowns (Kelkar, 2008; Kumar, 1987). Practically, gas network formulations can be coupled with large spectrum of friction factor calculation means for various flow regimes (Chen, 1979; Colebrook, 1939; GPSA, 2004; Menon, 2005; Moody, 1944; Weymouth, 1912). The generalized Newton-Raphson method is one of the most powerful solution schemes which can be implemented to solve many types of nonlinear systems of equations. Its significant drawback is its wellknown local (rather than global) convergence behaviordin which convergence is only guaranteed if the Newton Raphson scheme is initialized close-enough to the actual solution. This necessitates the availability of good initial guesses for all unknownsdwhich can become a daunting task for large and complex systems of equations
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
and unknowns. In order to overcome this obstacle, Ayala and Leong (2013) developed a robust and guess-free linear-pressure analog method to linearize and solve the system of equations written in terms of the p-formulation. The authors demonstrated that the linearpressure analog is able to reliably converge to a final solution without the need for any pressure or flow-rate initialization. More importantly, it was shown that the analog methodology could be utilized to aid the Newton-Raphson scheme by providing very good initializations (Ayala and Leong, 2013; Leong and Ayala, 2014). This study pursues to extend this previous work by presenting a linear-rate analog capable of effectively solve the q- (rather than the p-) formulation. The implementation of the q-formulation is practical and desirable because the resulting system of nodal mass balance equations are all linear when written in terms of flow ratedwhich constitute the large majority of network equations within the q-formulation. Other available alternatives for the solution of equations in terms of the q-formulation are the Hardy Cross (1930) and Linear Theory (Wood and Charles, 1972) methods. However, the implementation of the aforementioned solution techniques also requires the availability of good initial guesses. Some more recent studies have focused in speeding up the convergence of the classical methods such as Hardy Cross (e.g., Brkic, 2009) but little attention has been paid into eliminating the need for pressure or flow initial guesses from these protocols as presented in this study. In this work, the proposed linear rate analog linearizes the remaining non-linear loop equationsdleading to a matrix system whose coefficients are fully independent from flow or pressure initializations. Detailed case studies are presented to highlight the performance of the proposed method in solving horizontal inclined network cases. Numerical performance is compared against other linearization methodologies. It is shown that the proposed method represents a robust means to generate stable and fast solutions for the q-formulation which could be used to provide valuable guidance during gas pipe network design and optimization.
231
are flow rates at each of the pipe branches. The system of equations consists of a combination of nodal mass balance and loop equations. Nodal mass balance equations are the same as the ones presented in Equation (1), but written in terms of pipe flow rates. For an arbitrary node i, the nodal mass balance equation is expressed as follows:
X
±qij þ Si Di ¼ 0
(2)
j
where qij's are all the flows entering are leaving node i through pipe connections with neighboring j-nodes, and Si and Di are the external supplies and demands at node i, respectively. Flows entering the node are considered positive and flows leaving the node negative. Equation (2) is analogous to the first law of Kirchhoff for conservation of electrical charge in electrical circuits. In a pipe network of with N nodes, N-1 linear dependent equations can be generated via Equation (2). If “L” loops can be identified in the network system, additional loop equations can be established for each loop via the second law of Kirchhoff. Around every closed loop, the sum of all pressure drops around the geometrical enclosure must add up to zero: loop X ± p2i p2j ¼ 0
(3)
i;j
In this summation (Equation (3)), pressure drops that result in flow directions that align with prescribed loop direction are considered positive, and negative otherwise. According to the generalized gas pipe equation for the horizontal case (see Appendix A), one writes:
p2i p2j ¼
qij Cij
!2 ¼ Rij q2ij
(4)
2. Governing equations in gas network analysis
where Rij is the pipe resistance (Rij ¼ 1=Cij2 ). Substituting Equation (4) into Equation (3), loops equations could be alternatively recast as:
2.1. p-formulation
loop X
In the p-formulation, continuity equations for gas flow in conduits are written in terms of nodal pressures which are treated as primary unknowns. For a horizontal pipeline network, gas network equations are written at each node i considering all the connecting j-nodes via pipe branches as shown in Equation (1):
X
±Cij
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2i p2j þ Si Di ¼ 0
(1)
j
where Cij is the pipe conductivity, and Si and Di are flow sources or demands at node i, respectively. Flows entering the node are considered positive and flows leaving the node negative. Pipe conductivity calculations involves the calculation of the appropriate friction factor depending on the gas pipe equation of interest (see Appendix A). When the p-formulation is applied to a network of N nodes, only N-1 linearly independent equations can be obtained from Equation (1). Mathematical closure is achieved specifying pressure at one of the ‘N’ nodes in the system.
±Rij q2ij ¼ 0
(5)
ij
It should be reiterated at this point that the resulting system of equations generated using the q-formulation is mostly linear: “N-1” equations are already linear in flow rate (nodal mass conservation, Equation (2)) with additional “L” non-linear equations (loop equations, Equation (5)) because of the presence of squared flow (q2ij ) terms. In general, N >> L, and the system of equations would be fully linear if there were no closed loops in the network (L ¼ 0). 3. The linear-pressure analog formulation Ayala and Leong (2013) developed a linear-pressure analog scheme that linearized the p-form of the gas network equations. They showed that the non-linear pressure coefficient in the generalized gas flow equation (Equation (6)):
qij ¼ Cij
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2i p2j
(6)
can be rewritten as Equation (7): 2.2. q-formulation This work aims to solve governing equations for natural gas flow in conduits using the q-formulation. Primary unknowns in this case
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2i p2j ¼
sffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pi þ pj 2 pi pj ¼ 1 þ pi $ pi pj pi pj p 1 j
(7)
232
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Therefore, Equation (1) could be successfully linearized if written as Equation (8):
þ
X
Lij pi pj þ Si Di ¼ 0
tij ¼ (8)
j
where Lij is the linear-analog pipe conductivity defined as Equation (9):
Lij ¼ Cij
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1þ rij 1
(9)
where rij ¼ ppij is the pipe pressure ratio and i and j are the indices for the pipe upstream and downstream nodes, respectively. Equation (8) can be solved directly using any linear solver. While calculating nodal pressures iteratively, the pressure dependent term ri,j can be computed using results from last iteration. For the first iteration, Lij is made equal to any multiple of the actual pipe conductivity. Such multiple must be higher than one. The linearpressure analog demonstrated to be a robust and initial pressure and flow guess-free means of generating stable and reliable solutions for horizontal and inclined gas network systems (Ayala and Leong, 2013; Leong and Ayala, 2014).
(14)
where rij ¼ ppij . Since pi (upstream pipe pressure) is always larger than pj (downstream pressure), it follows that rij is always larger than 1 (pi > pj ). Thus, values of the tij analog transform variable must always result in real numbers larger than 0 and smaller than 1, as illustrated in Fig. 1. tij ¼ 0 corresponds to the limiting no-flow conditions at which pi ¼ pj (rij ¼ 1Þ and tij ¼ 1 represents the extreme condition pi >> pj (rij /∞Þ. For all practical purposes, tij is to be found within the approximate range 0 < tij < 0.5 because pipes seldom operate at rij's any larger than 2 given the excessive pressure (and thus economic) losses that would be involved. This is practical cut-off is also highlighted in Fig. 1. On the basis of Equation (12) and Equation (13), the linear-flow analog flow equation model proposed by this work then becomes:
pi pj ¼ lij qij
(15)
which readily allows to rewrite the pressure-drop closure equation around a loop (Kirchhoff's second law) in terms of linear pressure differences (Equation (16)) rather than squared pressure differences (Equation (3)): loop loop X X ± pi pj ¼ xij lij qij ¼ 0
4. A linear-rate analog formulation
i;j
The earlier successful development of linear-pressure analogs for the solution of gas networks motivated the exploration of linearization of the loop equations in q-formulation using a similar analog transformation. A linearized analog system could be expected to be more straightforward to solve because all nodal mass balance equations are already linear in terms of flow rate. Thus, linearization would only need to focus on loop equations (Equation (5)). Such linearization starts with Equation (4), which can be rewritten as Equation (10):
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q qffiffiffiffiffiffi ij ¼ Rij qij p2i p2j ¼ Cij
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 1 þ rij
(10)
In Equation (16), xij ¼ 1 if the pipe flow direction aligns with the prescribed positive direction of the loop, and xij ¼ 1 if it goes against it. At this stage, linearization of the loop equations has been accomplished. Once the linear system of equations is solved, nodal pressures and flow rates for each pipe branches can be readily obtained. Unlike Linear Theory, Hardy Cross, and even the Newton Raphson methods, the proposed solution not only linearizes the loop equations but also eliminates the dependencies of the resulting coefficient matrices on the primary unknowns (flow rates). More importantly, the major advantage of the proposed linear-rate analog is that no initial guesses of either nodal pressures nor flow rates are required to start the iterative scheme.
where i and j are the indices of upstream and downstream nodes, respectively. Thus:
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v ffi u u pi þ pj q u pi pj ¼ ij t Cij pi pj
(11)
Rearranging Equation (11), one may obtain Equation (12):
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 3 u " sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi # u pi pj 1 1 2 u 7 6 t 5qij ¼ pi pj ¼ 4 1 q Cij C 1 þ rij ij ij p þp 2
i
(12)
j
from where, a ‘linear-rate analog resistance’ (lij Þ can be defined as shown in Equation (13):
lij ¼
1 Cij
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 1 ¼ t 1 þ rij Cij ij
(13)
where tij, the ‘analog linear-resistance transform’, is being defined as:
(16)
ij
Fig. 1. Behavior of the analog linear-resistance transform tij versus rij.
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
233
4.1. Proposed solution strategy Figure 2 illustrates the successive substitution scheme recommended for the proposed linear analog solution processes. The protocol is initialized by setting lij ¼ C1ij in all loop equationsdi.e, with initial tij equal to 1 (to ¼ 1Þ for all pipes during the first iteration only. It should be noted that this initialization rescinds any need to know initial flow rates. The (non-)impact of selecting a different initial tij value (other than tij ¼ to ¼ 1) is discussed in a later section in this paper. The linearized set can be straightforwardly solved via any linear solver of simultaneous equations, and flow rates can be calculated as a result. With flow rates and the specified nodal pressure in the problem, nodal pressure distribution can be directly calculated along with pipe pressure ratios (rij ). Once rij ’s become available, linear-rate analog resistivities (lij ) are updated via Equation (13). Iteration continues until the desired flow rate tolerance is achieved. It should be kept in mind that the sign of xij should always be updated if a pipe flow reversal occurs during iterations (negative flow rates). 5. Case studies 5.1. Case study 1: horizontal pipe network system In this section, the linear rate analog is implemented to solve the same horizontal pipeline network presented by Ayala and Leong (2013) of nine nodes and twelve branches. As displayed in Fig. 3, four loops can be identified in the network and the clockwise direction is prescribed as positive. Pressure at Node 9 is anchored (specified) at 130 psia. Table 1 summarizes loops which are specified using the passing nodes along the positive direction. It is noted that the discussions of this paper are completely general, but these case studies utilize AGA fully-turbulent friction factor (AGA, 1965) calculation for illustration purposes. Additional network details are listed in Appendix B.
Fig. 2. Proposed workflow for the linear-rate analog algorithm.
Demand 2 MMSCFD
2
Pipe#1 (1,2)
Loop #1
Demand 2 MMSCFD
Demand 1 MMSCFD
Loop #3 Pipe#11 (7,8)
Pipe#5 (3,6) Pipe#8 (5,6)
Loop #4
8
Pipe#12 (8,9)
Demand 2 MMSCFD Fig. 3. Network topology and loops of Case 1.
6
Demand 2 MMSCFD
Pipe#10 (6,9)
Pipe#6 (4,5)
7
Demand 2 MMSCFD
3
Loop #2
5 Pipe#9 (5,8)
4 Pipe#7 (4,7)
Demand 3 MMSCFD
Pipe#3 (2,3)
Pipe#4 (2,5)
1 Pipe#2 (1,4)
Supply 16 MMSCFD
9
Demand 2 MMSCFD psf = 130 psia
234
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Table 1 Summary of the loops in Case 1.
Linearized loop equations
Nodal mass balance equations
Fig. 4 shows the resulting system of equations generated by the linear-rate analog method in a matrix form. Twelve linear equations are solved to update twelve flow rates in the pipeline branches. The first eight equations are nodal mass balance equations while the other four equations are linearized loop equations. This system of equations can be solved utilizing any linear solver. For comparison purposes, three linear solvers are used to solve this problem: the linear-rate analog (this work), linear-pressure analog (Ayala and Leong, 2013), and the Linear theory method (Wood and Charles, 1972). For the linear theory solution, a random distributed gas flow rates between 1 and 2 MMSCF/day is used as initial flow guess because it is required for linear-theory iterations to begin. Table 2 and Table 3 show that the prediction of node pressures and flow rates are identical by implementing any of the three aforementioned methods. Fig. 5 shows that it takes 38 and 34 iterations and 6.705 ms and 6.399 ms of CPU time for the linear theory and linear-pressure analog, respectively, to achieve the prescribed flow rate tolerance (1 104 SCF/day). Linear-flow analog takes, in turn, only 19 iterations and 4.836 ms of CPU time to achieve the same convergence of 1 104 SCF/day. The linearrate analog also exhibits a more robust performance than linear theory for this case study. An important reason for this is that the resulting coefficient matrix is not explicitly dependent of flow rates in the linear-rate case, which significantly speeds up convergence and eliminates the need for flow-rate initializations. In addition, the linear-rate analog takes fewer iterations than even its linearpressure counterpart, for the same identical case study. This can be explained by comparing the non-linear complexity of both formulations: the system of equations generated by the q-formulation contains less nonlinear equations than that of p-formulation.
Table 2 Summary of node pressure prediction using different methods, Case Study 1. Nodal pressure, psia Node
Linear-rate analog
Linear-pressure analog (Ayala and Leong, 2013)
Linear theory
1 2 3 4 5 6 7 8 9
536.112 378.833 328.874 351.089 264.165 198.796 300.382 187.528 130
536.112 378.833 328.874 351.089 264.165 198.796 300.382 187.528 130
536.112 378.833 328.874 351.089 264.165 198.796 300.382 187.528 130
Table 3 Summary of flow rate predictions using different methods, Case Study 1. Gas rate, SCF/day Inlet
Outlet
Linear-rate analog
Linear-pressure analog (Ayala and Leong, 2013)
Linear theory
1 1 2 2 3 4 4 5 5 6 7 8
2 4 3 5 6 5 7 6 8 9 8 9
7 736 768.630 8 263 231.370 3 834 954.289 1 901 814.340 1 834 954.289 1 619 725.416 3 643 505.954 1 218 428.476 1 303 111.281 1 053 382.765 1 643 505.954 946 617.235
7 736 768.630 8 263 231.370 3 834 954.289 1 901 814.340 1 834 954.289 1 619 725.416 3 643 505.954 1 218 428.476 1 303 111.281 1 053 382.765 1 643 505.954 946 617.235
7 736 768.630 8 263 231.370 3 834 954.289 1 901 814.340 1 834 954.289 1 619 725.416 3 643 505.954 1 218 428.476 1 303 111.281 1 053 382.765 1 643 505.954 946 617.235
-1
-1
0
0
0
0
0
0
0
0
0
0
q1
-16000000
1
0
-1
-1
0
0
0
0
0
0
0
0
q2
2000000
0
0
1
0
-1
0
0
0
0
0
0
0
q3
2000000
0
1
0
0
0
-1
-1
0
0
0
0
0
q4
3000000
0
0
0
1
0
1
0
-1
-1
0
0
0
q5
1000000
0
0
0
0
1
0
0
1
0
-1
0
0
0
0
0
0
0
0
1
0
0
0
-1
0
q7
2000000
0
0
0
0
0
0
0
0
1
0
1
-1
q8
2000000
-λ1,4
0
λ2,5
0
- λ4,5
0
0
0
0
0
0
q9
0
0
0
λ2,3
- λ2,5
λ3,6
0
0
- λ5,6
0
0
0
0
q10
0
0
0
0
0
0
λ4,5
-λ4.7
0
λ5,8
0
- λ7,8
0
q11
0
0
0
0
0
0
0
0
λ5,6
- λ5,8
λ6,9
0
- λ8,9
q12
0
λ1,2
Fig. 4. Coefficient matrix of linear rate analog of Case Study 1.
×
q6
=
2000000
Fig. 5. Comparison of numerical performance of various linear solution techniques.
Fig. 6. Performance behavior - the analog linear-resistance transform t of the pipe branches.
236
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Fig. 7. Performance behavior - Gas flow rate solution.
Fig. 6, Fig. 7 and Fig. 8 display the convergence performance behavior of the linear-rate analog in terms of the evolution of the values of tij , flow rate, and pressure, respectively. As indicated in Fig. 2, the first iteration starts with tij ¼ to ¼ 1 for all pipes. It is noted in all these figures that convergence behavior of the proposed method is very smooth for all variables involved. Values of analog linear-resistance transform (tij ), flow rates, and pressure distributions become quite stable after the sixth iteration and updates for pipe branch flow rates and node pressures become small after few iterations. Additional numerical experiments were carried to quantify the impact of using different initial tij values (t0) on convergence behavior (other than the recommended tij ¼ 1 ¼ to used for the first iteration). As shown in Fig. 9, dramatically different initial t0 values (t0 < 0.01) were utilized compared to the to ¼ 1 choice. It is noted that convergence behavior of tij (and thus lij , since Cij is mostly or relatively constant depending on the choice of gas equation) is completely unaffected by the choice of to for the first iteration. This is due to the structure of the loop equations. For example, the expanded loop equation for Loop 1 (see Fig. 3) is shown below as Equation (17):
Equation (13), Equation (17) can be written as Equation (18):
t1;2
C1;2
q1
t1;4 C1;4
q2 þ
t2;5 C2;5
q4
t4;5 C4;5
q6
¼0
(18)
Therefore, for the very first iteration, this equation becomes:
t0
1 1 1 1 q q þ q q C1;2 1 C1;4 2 C2;5 4 C4;5 6
¼0
(19)
which demonstrates that using any value of t0 would not alter results from the first iteration given that such value cancels out of the equation (Equation (19)). Hence, unlike the linear-pressure analog where different initial linear-conductivity initializations can be shown to have some impact on convergence, it is demonstrated that convergence behavior of the proposed linear-rate analog is not sensitive to the t0 values used in the first iteration, and that the recommended use of lij ¼ 1=Cij (to ¼ 1) would suffice.
5.2. Case study 2: inclined pipe network system
l12 q1 l14 q2 þ l25 q4 l45 q6 ¼ 0
(17)
Expressing the linear rate analog resistivity terms using
In order to consider the elevation of the pipeline network, Equation (4) is modified written as Equation (20) (Ferguson, 1936):
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
237
Fig. 8. Performance behavior - Nodal pressure solution.
p2i esij p2j ¼
qij Cij
the inclined pipe Equation (20) can be rewritten for the case of a linear-q formulation as follows:
!2 (20)
where sij is a term which characterizes the elevation different between the upstream and downstream nodes, which is defined in Equation (21) (Ayala, 2013):
sij ¼ 0:0375
gg D H zave TAve
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 sij u " sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi # u pi e 2 pj 1 2 7 61 u q pi e pj ¼ 4 t ¼ 1 q 5 ij sij Cij C 1 þ rij ij ij 2 pi þ e pj 2
sij 2
¼
1 t Cij ij
(21)
(24)
In this case, the pipeline conductivity is written in terms of equivalent pipe lengths, as shown in Equation (22) (Ayala, 2013):
where the linear-rate analog resistance (lij Þ and the analog linearresistance transform (tij) are again defined for the inclined case as:
Ef sG pTscsc 1 d2:5 Cij ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi 0:5 gg Tave zave fF Le
1 lij ¼ Cij
(22)
In Equation (22), Le is the equivalent pipeline length, in ft, which can be calculated as (Ayala, 2013):
Le ¼
ðesij 1Þ L sij
(23)
Following the development presented in Equations (10)e(14),
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 1 ¼ t 1 þ rij Cij ij
(25)
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 tij ¼ 1 1 þ rij
(26) sij
but with rij ¼ spij i . Note the presence of the ‘e 2 ’ in denominator of 2 pj the rij ratio foreinclined pipes. Therefore, the linear-flow analog flow equation model proposed by this work for the inclined case is:
238
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Fig. 9. Impact of initial values transformation variable t on the iteration efficacy.
sij
pi e 2 pj ¼ lij qij
(27)
Or,
sij pi pj ¼ lij qij þ e 2 1 pj
(28)
As a result, Kirchhoff's second law (Equation (16)) is rewritten as follows for the inclined case as follows: loop loop X h sij i X ± pi pj ¼ xij lij qij þ e 2 1 pj ¼ 0 i;j
(29)
ij
or, loop X ij
xij lij qij ¼
loop X
sij
xij 1 e 2 pj
(30)
ij
Again, xij equals to positive one if the flow direction agrees with the prescribed loop positive direction, or equals to negative one if it goes against. The strategy is similar to the horizontal case, except that the value of the right-hand side (RHS) of Equation (30) needs to be updated at every iteration other than the first. For the first
iteration, tij ¼ 1 ¼ to and the RHS of Equation (30) is ignored (as if it were a horizontal system) to startup the iteration process. Linear-rate analog method is now implemented to solve the inclined pipeline network proposed by Ayala and Leong (2013), and displayed in Fig. 10. This network system is similar to Case study 1 (four loops and nine nodes) but nodes at located at different elevations. Results obtained via the implementation of the linear-rate analog scheme are presented in Table 4 and Table 5 in terms of nodal pressure and pipeline branch flow rates, respectively. It is noteworthy that the solution of linear rate analog shows an exact agreement with the reported results generated by Ayala and Leong using linear pressure analog. As shown in Fig. 11, however, the linear-rate analog takes 19 iterations to achieved the prescribed flow rate tolerance of 1 104 SCF/day, while linear-pressure analog takes 34 iterations and linear theory takes 39 iterations to achieve the same level of tolerance. Also, the required CPU time for linear-rate analog is smaller compared to the other linearization techniques. It should be added that the linear theory method was not able to originally converge for this case study until all known flow directions were explicitly specified (which was done in order to generate the comparison in Fig. 11). Neither linear-pressure nor linear-rate analogs require for any branch flow directions to be known to ensure convergence.
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
239
Demand 2 MMSCFD
7
Δh (7,8) = 600 ft
Pipe#11 (7,8)
Pipe#8 (5,6)
8
Loop #4
Δh (8,9) = -1000 ft
Pipe#12 (8,9)
Pipe#5 (3,6)
Δh (3,6) = - 800 ft
Δh (2,5) = - 400 ft
Pipe#4 (2,5)
Δh (5,6) = 1200 ft
6
Demand 2 MMSCFD
Pipe#10 (6,9)
Loop #3
Loop #2
Δh (6,9) = -1600 ft
Demand 1 MMSCFD
Demand 2 MMSCFD
3
Pipe#3 (2,3)
5
Pipe#6 (4,5) Pipe#7 (4,7)
Δh (4,7) = - 800 ft
Pipe#1 (1,2)
Loop #1
Δh (2,3) = 1200 ft
2
Δh (4,5) = -800 ft
4
Demand 2 MMSCFD
Δh (1,2) = -1200 ft
Pipe#9 (5,8)
Demand 3 MMSCFD
Pipe#2 (1,4)
Δh (1,4) = - 800 ft
1
Δh (5,8) = 600 ft
Supply 16 MMSCFD
9
Demand 2 MMSCFD psf = 130 psia
Demand 2 MMSCFD Fig. 10. Network topology and identified loops for inclined pipeline network system.
6. Concluding remarks
Table 4 Summary of the node pressure prediction, Case study 2. Node pressure, psia Node
Linear-rate analog
Linear-pressure analog (Ayala and Leong, 2013)
Linear theory
1 2 3 4 5 6 7 8 9
524.759 380.602 319.392 346.358 266.612 191.583 302.400 185.047 130
524.759 380.602 319.392 346.358 266.612 191.583 302.400 185.047 130
524.759 380.602 319.392 346.358 266.612 191.583 302.400 185.047 130
Table 5 Summary of the gas flow rate prediction, Case study 2. Gas rate, SCF/day Inlet
outlet
Linear-rate analog
Linear-pressure analog (Ayala and Leong, 2013)
Linear theory
1 1 2 2 3 4 4 5 5 6 7 8
2 4 3 5 6 5 7 6 8 9 8 9
7 8 3 1 1 1 3 1 1 1 1 9
7 8 3 1 1 1 3 1 1 1 1 9
7 8 3 1 1 1 3 1 1 1 1 9
742 346.142 257 653.858 810 225.731 932 120.412 810 225.731 612 367.344 645 286.513 231 468.278 313 019.478 041 694.009 645 286.513 58 305.991
742 346.142 257 653.858 810 225.731 932 120.412 810 225.731 612 367.344 645 286.513 231 468.278 313 019.478 041 694.009 645 286.513 58 305.991
742 346.142 257 653.858 810 225.731 932 120.412 810 225.731 612 367.344 645 286.513 231 468.278 313 019.478 041 694.009 645 286.513 58 305.991
In this work, a linear-rate analog formulation is proposed to solve the natural gas transportation in pipeline networks problems using a q- (flow-) formulation. Linear-rate analog effectively linearizes the loop equations formulated by Kirchhoff's second law of a gas pipeline network system via the use of linear-rate analog pipeline resistivities. Through case studies, the linear-rate analog is shown to be a robust, stable and initial guess-free solution technique for natural gas pipeline network analysis. The following additional observations may be drawn: Linear-rate analog can be successfully implemented to solve both horizontal and inclined pipeline network cases. The solutions of linear-rate analog agree well with known reported results in terms of node pressures and branch flow rates for both cases. The numerical performance of the linear-rate analog method has been compared against two other avalaible linearization algorithms, linear theory and linear-pressure analog method, by solving the same pipeline network cases. It is shown that the linear-rate analog requires fewer iteration for convergence and smaller CPU time than the other two solution techniques for all cases under consideration. The linear-rate analog method is an initial-guess free technique (no pressure or flow rate guesses are needed to start iterations) that linearizes the loop equations and constructs a coefficient matrix that is independent of primary unknowns (flow rate). Computational savings associated with the use of linear-rate analog method is expected to grow considerably as the size of the network system increases because most of the equations deployed by the q-formulation are linear equations.
240
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Fig. 11. Comparison of numerical performances of linear-rate analog against linear theory and linear-pressure analog e Case Study 2.
Nomenclature AGA Cij Di d e fF DH L pi qij r rG R Rij s Si T z
American Gas Association pipeline conductivity of a pipe branch connected by nodes i and j SCF/day-psi gas demand at an arbitrary node i SCF/day pipeline diameter inch pipeline toughness ft Fanning friction factor dimensionless elevation difference between nodes ft pipeline length ft pressure of an arbitrary node i psia gas flow rate of a pipe branch connected by nodes i and j SCF/day pipeline node pressure ratio dimensionless pipeline specific resistivity psia2-inch5/SCFD2-ft gas constant, with a magnitude of 10.732 psia-ft3lbmole1- R1 pipeline resistivity of pipe branch connected by node i and j psia2/SCF2/day2 pipe elevation parameter dimensionless gas supply at an arbitrary node i SCF/day temperature R gas compressibility factor dimensionless
Appendix A. Summary of gas flow equations
Table A-1 Summary of the gas equations (Modified from Ayala, 2013). Generalized gas equation in conduits: qij ¼ Cij Ef sG
Tsc psc
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2i esij p2j ,
Cij ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p1ffiffiffi gg Tave zave
d2:5 0:5 fF Le
Gas equation General equation
Definition of friction factor ! p1ffiffiffi ¼ 4:0 log10 fF
Weymouth
fF ¼ 0:008 d1=3
Panhandle-A
fF ¼
0:01923
fF ¼
0:00359
qgSC gg d
Panhandle-B
qgSC gg d
AGA Fully Turbulent
3:7 þ
0:03922
3:7d e
where Re is the Reynolds number, qgSC is in SCF/day, L is in ft or mile, d is in inch, p is in psia and T is in R. The unit conversion factor sG varies based on the unit of the pipe length. For L in ft, sG ¼ 2818, sG,w ¼ 31 508, sG,PA ¼ 44 400, sG,PA ¼ 58 328; For L in miles, sG ¼ 3 ¼ G,w ¼ 4 ¼ sG,PA ¼ 435.98, sG,PA ¼ 736.
Greek
g l s x t
gas specific gravity dimensionless Linear-rate analog pipeline resistivity psi/SCF/day gas equation unit conversion factor directional sign term of a close pipeline loop transformation variable
Subscripts ave average value e equivalent i branch upstream node j branch downstream node sc standard condition 0 first iteration
5:02 pffiffiffi Re fF
0:1461
p1ffiffiffi ¼ 4 log10 fF
e d
Appendix B. Network details e Case Studies
Table B-1 Gas properties and general data (Case Study 1 and 2). Base pressure (SC)
14.7 psia
Base temperature (SC) Average gas compressibility Average gas gravity T average Pipe roughness Pipe efficiency
60.0 F 0.9 0.69 535 R 0.0018 in 1
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
241
Table B-2 Pipe and supply/demand data (Case Study 1 and 2).
References American Gas Association AGA, 1965. Steady Flow in Gas Pipelines: Testing,/ Measuremen,/Behavior, and Computation. Institute of Gas Technology, New York, NY, U.S. Ayala, H.,L.F., 2013. Chapter 22: Transportation of crude oil, natural gas, and petroleum products. In: Petroleum Refining and Natural Gas Processing. ASTM International, West Conshohocken, PA, USA, pp. 549e575. Ayala, L.F., Leong, C.Y., 2013. A robust linear-pressure analog for the analysis of a natural gas transportation networks. J. Nat. Gas Sci. Eng. 174e184. Brki c, D., 2009. An improvement of Hardy Cross method applied on looped spatial. Appl. Energy 86, 1290e1300. Chen, N.H., 1979. An explicit equation for friction factor in pipe. Ind. Eng. Chem. Fund. 18, 296e297. Colebrook, C.F., 1939. Turbulent flow in pipes, with particular reference to the transition region between smooth and rough pipe laws. J. Inst. Civ. Eng. 12, 133e156. Cross, H., 1930. Analysis of continuous frames by distributing fixed-end moments. Proc. Am. Soc. Civ. Eng. 56 (5), 919e928. Energy Information Administration EIA, 2013. Technically Recoverable Shale Oil and Shale Gas Resources: an Assessment of 137 Shale Formations in 41 Countries
outside the United States. EIA Independent Statistics and Analysis. U.S. Energy Information Administration, Washington DC. U.S. Ferguson, J.W., 1936. A mathematical analysis of the effect of differences on flow formulae for gas pipe lines. Gas. Age 78, 72e76. Menon, E.S., 2005. Gas Pipeline Hydraulics. Taylor and Francis Group, Boca Raton, FL, U.S. Gas Processors Supplier Association GPSA, 2004. “Fluid Flow and Piping.” GPSA Engineering Data Book. 12. Gas Processors Supplier Association, Tulsa, OK, U.S. Kelkar, M., 2008. Natural Gas Production Engineering. Penn Well Corporation, Tulsa, OK, U.S. Kumar, S., 1987. Gas production Engineering. Gulf Pub. Co., Houston. Leong, C.Y., Ayala, L.F., 2014. Hybrid approach using linear analogs for gas network simulation with multiple components. Oil Gas Facil. J. 3, 76e88. Makogon, Y.F., Holditch, S.A., Makogon, T.Y., 2007. Natural gas-hydrates-a potential energy source. J. Pet. Sci. Eng. 56, 14e31. Moody, L.F., 1944. Friction factors for pipe flow. Trans. ASME 66, 671e684. Weymouth, T.R., 1912. Problems in natural gas engineering. Trans. ASME 34, 185e206. Wood, D.A., Charles, C.O., 1972. Hydraulic network analysis using linear theory. ASCE J. Hydraul. Div. 1157e1170. Wood, D.A., Mokhatab, S., Economides, M.J., 2008. Technology options for securing markets for remote gas. In: Proceedings of 87th Annual Convention of the Gas Processors Association. GPA, Grapevine, Texas, US.
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
A linear-rate analog approach for the analysis of natural gas transportation networks Qian Sun, Luis F. Ayala* John and Willie Leone Family Department of Energy and Mineral Engineering, The Pennsylvania State University, University Park, PA 16802, United States
a r t i c l e i n f o
a b s t r a c t
Article history: Received 26 October 2016 Received in revised form 18 January 2017 Accepted 27 March 2017 Available online 6 April 2017
Modeling natural gas transportation networks poses a number of challenges due to the significant nonlinearities associated to the governing equations. Pressure (nodal or p-) formulations and flow (nodalloop or q-) formulations are the most commonly deployed approaches used to formulate gas flow network models. They treat nodal pressures and pipe branch flow rates as primary unknowns, respectively, and the Newton-Raphson method is the typical choice used to solve the resulting system of pipe network equations. The major disadvantage of Newton-Raphson methods is their tendency to hopelessly diverge when good initializations for the unknown pressure and flow variables are not available. In order to overcome this drawback, a linear-pressure analog approach, applicable to the p-formulation, was recently proposed to formulate an initial-guess-free solution protocol. However, solving the q-formulation rather than the p-formulation would be an alternate, practical, and desirable way to study gas flow in pipeline networks because the system of equations is predominantly lineardwith the exception of the loop equations. Linear Theory and Hardy Cross methods have been used in the past solve such qformulation. Unfortunately, these methods are not initial-guess-free and have been shown to become potentially unstable and thus inefficient because their convergence is also strongly associated with the availability of good initial guesses. This work proposes a linear-rate analog method capable of effectively solving the q-formulation. Through case studies, the proposed linear-rate analog method is shown to be a robust and initial-guess-free solution scheme to effectively model horizontal and inclined natural gas pipeline networks. © 2017 Elsevier B.V. All rights reserved.
Keywords: q-formulation Gas flow in conduits Networks Linear-rate analog
1. Introduction Reliable modeling of natural gas transportation networks continues to increase in importance as a result of increased energy demands and production from large unconventional gas fields (EIA, 2013; Makogon et al., 2007; Wood et al., 2008). Solving a natural gas pipeline network system is challenging due to the high nonlinearity of the governing equations of gas flow in conduits. Therefore, the development of robust mathematical methods to model gas flow in pipeline becomes of utmost importance for the study of optimized design and performance. Governing equations of gas pipe flow in networks are derived from energy conservation statements for steady state conditions (Ayala, 2013). Typically, pipe networks consist of nodes and pipe branches. External demands and supplies can be present at the
* Corresponding author. E-mail address: [email protected] (L.F. Ayala). http://dx.doi.org/10.1016/j.jngse.2017.03.027 1875-5100/© 2017 Elsevier B.V. All rights reserved.
nodes. In these networks, closed circuits formed by pipe branches are known as pipe loops. The formulation of gas flow network equations can be categorized as pressure (p-), pipe flow (q-) or loop-rate correction (Dq-) approaches, where node pressures, branch flow rates or corrections to flow rates within a loop, respectively, are treated as primary unknowns (Kelkar, 2008; Kumar, 1987). Practically, gas network formulations can be coupled with large spectrum of friction factor calculation means for various flow regimes (Chen, 1979; Colebrook, 1939; GPSA, 2004; Menon, 2005; Moody, 1944; Weymouth, 1912). The generalized Newton-Raphson method is one of the most powerful solution schemes which can be implemented to solve many types of nonlinear systems of equations. Its significant drawback is its wellknown local (rather than global) convergence behaviordin which convergence is only guaranteed if the Newton Raphson scheme is initialized close-enough to the actual solution. This necessitates the availability of good initial guesses for all unknownsdwhich can become a daunting task for large and complex systems of equations
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
and unknowns. In order to overcome this obstacle, Ayala and Leong (2013) developed a robust and guess-free linear-pressure analog method to linearize and solve the system of equations written in terms of the p-formulation. The authors demonstrated that the linearpressure analog is able to reliably converge to a final solution without the need for any pressure or flow-rate initialization. More importantly, it was shown that the analog methodology could be utilized to aid the Newton-Raphson scheme by providing very good initializations (Ayala and Leong, 2013; Leong and Ayala, 2014). This study pursues to extend this previous work by presenting a linear-rate analog capable of effectively solve the q- (rather than the p-) formulation. The implementation of the q-formulation is practical and desirable because the resulting system of nodal mass balance equations are all linear when written in terms of flow ratedwhich constitute the large majority of network equations within the q-formulation. Other available alternatives for the solution of equations in terms of the q-formulation are the Hardy Cross (1930) and Linear Theory (Wood and Charles, 1972) methods. However, the implementation of the aforementioned solution techniques also requires the availability of good initial guesses. Some more recent studies have focused in speeding up the convergence of the classical methods such as Hardy Cross (e.g., Brkic, 2009) but little attention has been paid into eliminating the need for pressure or flow initial guesses from these protocols as presented in this study. In this work, the proposed linear rate analog linearizes the remaining non-linear loop equationsdleading to a matrix system whose coefficients are fully independent from flow or pressure initializations. Detailed case studies are presented to highlight the performance of the proposed method in solving horizontal inclined network cases. Numerical performance is compared against other linearization methodologies. It is shown that the proposed method represents a robust means to generate stable and fast solutions for the q-formulation which could be used to provide valuable guidance during gas pipe network design and optimization.
231
are flow rates at each of the pipe branches. The system of equations consists of a combination of nodal mass balance and loop equations. Nodal mass balance equations are the same as the ones presented in Equation (1), but written in terms of pipe flow rates. For an arbitrary node i, the nodal mass balance equation is expressed as follows:
X
±qij þ Si Di ¼ 0
(2)
j
where qij's are all the flows entering are leaving node i through pipe connections with neighboring j-nodes, and Si and Di are the external supplies and demands at node i, respectively. Flows entering the node are considered positive and flows leaving the node negative. Equation (2) is analogous to the first law of Kirchhoff for conservation of electrical charge in electrical circuits. In a pipe network of with N nodes, N-1 linear dependent equations can be generated via Equation (2). If “L” loops can be identified in the network system, additional loop equations can be established for each loop via the second law of Kirchhoff. Around every closed loop, the sum of all pressure drops around the geometrical enclosure must add up to zero: loop X ± p2i p2j ¼ 0
(3)
i;j
In this summation (Equation (3)), pressure drops that result in flow directions that align with prescribed loop direction are considered positive, and negative otherwise. According to the generalized gas pipe equation for the horizontal case (see Appendix A), one writes:
p2i p2j ¼
qij Cij
!2 ¼ Rij q2ij
(4)
2. Governing equations in gas network analysis
where Rij is the pipe resistance (Rij ¼ 1=Cij2 ). Substituting Equation (4) into Equation (3), loops equations could be alternatively recast as:
2.1. p-formulation
loop X
In the p-formulation, continuity equations for gas flow in conduits are written in terms of nodal pressures which are treated as primary unknowns. For a horizontal pipeline network, gas network equations are written at each node i considering all the connecting j-nodes via pipe branches as shown in Equation (1):
X
±Cij
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2i p2j þ Si Di ¼ 0
(1)
j
where Cij is the pipe conductivity, and Si and Di are flow sources or demands at node i, respectively. Flows entering the node are considered positive and flows leaving the node negative. Pipe conductivity calculations involves the calculation of the appropriate friction factor depending on the gas pipe equation of interest (see Appendix A). When the p-formulation is applied to a network of N nodes, only N-1 linearly independent equations can be obtained from Equation (1). Mathematical closure is achieved specifying pressure at one of the ‘N’ nodes in the system.
±Rij q2ij ¼ 0
(5)
ij
It should be reiterated at this point that the resulting system of equations generated using the q-formulation is mostly linear: “N-1” equations are already linear in flow rate (nodal mass conservation, Equation (2)) with additional “L” non-linear equations (loop equations, Equation (5)) because of the presence of squared flow (q2ij ) terms. In general, N >> L, and the system of equations would be fully linear if there were no closed loops in the network (L ¼ 0). 3. The linear-pressure analog formulation Ayala and Leong (2013) developed a linear-pressure analog scheme that linearized the p-form of the gas network equations. They showed that the non-linear pressure coefficient in the generalized gas flow equation (Equation (6)):
qij ¼ Cij
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2i p2j
(6)
can be rewritten as Equation (7): 2.2. q-formulation This work aims to solve governing equations for natural gas flow in conduits using the q-formulation. Primary unknowns in this case
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2i p2j ¼
sffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pi þ pj 2 pi pj ¼ 1 þ pi $ pi pj pi pj p 1 j
(7)
232
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Therefore, Equation (1) could be successfully linearized if written as Equation (8):
þ
X
Lij pi pj þ Si Di ¼ 0
tij ¼ (8)
j
where Lij is the linear-analog pipe conductivity defined as Equation (9):
Lij ¼ Cij
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1þ rij 1
(9)
where rij ¼ ppij is the pipe pressure ratio and i and j are the indices for the pipe upstream and downstream nodes, respectively. Equation (8) can be solved directly using any linear solver. While calculating nodal pressures iteratively, the pressure dependent term ri,j can be computed using results from last iteration. For the first iteration, Lij is made equal to any multiple of the actual pipe conductivity. Such multiple must be higher than one. The linearpressure analog demonstrated to be a robust and initial pressure and flow guess-free means of generating stable and reliable solutions for horizontal and inclined gas network systems (Ayala and Leong, 2013; Leong and Ayala, 2014).
(14)
where rij ¼ ppij . Since pi (upstream pipe pressure) is always larger than pj (downstream pressure), it follows that rij is always larger than 1 (pi > pj ). Thus, values of the tij analog transform variable must always result in real numbers larger than 0 and smaller than 1, as illustrated in Fig. 1. tij ¼ 0 corresponds to the limiting no-flow conditions at which pi ¼ pj (rij ¼ 1Þ and tij ¼ 1 represents the extreme condition pi >> pj (rij /∞Þ. For all practical purposes, tij is to be found within the approximate range 0 < tij < 0.5 because pipes seldom operate at rij's any larger than 2 given the excessive pressure (and thus economic) losses that would be involved. This is practical cut-off is also highlighted in Fig. 1. On the basis of Equation (12) and Equation (13), the linear-flow analog flow equation model proposed by this work then becomes:
pi pj ¼ lij qij
(15)
which readily allows to rewrite the pressure-drop closure equation around a loop (Kirchhoff's second law) in terms of linear pressure differences (Equation (16)) rather than squared pressure differences (Equation (3)): loop loop X X ± pi pj ¼ xij lij qij ¼ 0
4. A linear-rate analog formulation
i;j
The earlier successful development of linear-pressure analogs for the solution of gas networks motivated the exploration of linearization of the loop equations in q-formulation using a similar analog transformation. A linearized analog system could be expected to be more straightforward to solve because all nodal mass balance equations are already linear in terms of flow rate. Thus, linearization would only need to focus on loop equations (Equation (5)). Such linearization starts with Equation (4), which can be rewritten as Equation (10):
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q qffiffiffiffiffiffi ij ¼ Rij qij p2i p2j ¼ Cij
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 1 þ rij
(10)
In Equation (16), xij ¼ 1 if the pipe flow direction aligns with the prescribed positive direction of the loop, and xij ¼ 1 if it goes against it. At this stage, linearization of the loop equations has been accomplished. Once the linear system of equations is solved, nodal pressures and flow rates for each pipe branches can be readily obtained. Unlike Linear Theory, Hardy Cross, and even the Newton Raphson methods, the proposed solution not only linearizes the loop equations but also eliminates the dependencies of the resulting coefficient matrices on the primary unknowns (flow rates). More importantly, the major advantage of the proposed linear-rate analog is that no initial guesses of either nodal pressures nor flow rates are required to start the iterative scheme.
where i and j are the indices of upstream and downstream nodes, respectively. Thus:
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v ffi u u pi þ pj q u pi pj ¼ ij t Cij pi pj
(11)
Rearranging Equation (11), one may obtain Equation (12):
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 3 u " sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi # u pi pj 1 1 2 u 7 6 t 5qij ¼ pi pj ¼ 4 1 q Cij C 1 þ rij ij ij p þp 2
i
(12)
j
from where, a ‘linear-rate analog resistance’ (lij Þ can be defined as shown in Equation (13):
lij ¼
1 Cij
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 1 ¼ t 1 þ rij Cij ij
(13)
where tij, the ‘analog linear-resistance transform’, is being defined as:
(16)
ij
Fig. 1. Behavior of the analog linear-resistance transform tij versus rij.
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
233
4.1. Proposed solution strategy Figure 2 illustrates the successive substitution scheme recommended for the proposed linear analog solution processes. The protocol is initialized by setting lij ¼ C1ij in all loop equationsdi.e, with initial tij equal to 1 (to ¼ 1Þ for all pipes during the first iteration only. It should be noted that this initialization rescinds any need to know initial flow rates. The (non-)impact of selecting a different initial tij value (other than tij ¼ to ¼ 1) is discussed in a later section in this paper. The linearized set can be straightforwardly solved via any linear solver of simultaneous equations, and flow rates can be calculated as a result. With flow rates and the specified nodal pressure in the problem, nodal pressure distribution can be directly calculated along with pipe pressure ratios (rij ). Once rij ’s become available, linear-rate analog resistivities (lij ) are updated via Equation (13). Iteration continues until the desired flow rate tolerance is achieved. It should be kept in mind that the sign of xij should always be updated if a pipe flow reversal occurs during iterations (negative flow rates). 5. Case studies 5.1. Case study 1: horizontal pipe network system In this section, the linear rate analog is implemented to solve the same horizontal pipeline network presented by Ayala and Leong (2013) of nine nodes and twelve branches. As displayed in Fig. 3, four loops can be identified in the network and the clockwise direction is prescribed as positive. Pressure at Node 9 is anchored (specified) at 130 psia. Table 1 summarizes loops which are specified using the passing nodes along the positive direction. It is noted that the discussions of this paper are completely general, but these case studies utilize AGA fully-turbulent friction factor (AGA, 1965) calculation for illustration purposes. Additional network details are listed in Appendix B.
Fig. 2. Proposed workflow for the linear-rate analog algorithm.
Demand 2 MMSCFD
2
Pipe#1 (1,2)
Loop #1
Demand 2 MMSCFD
Demand 1 MMSCFD
Loop #3 Pipe#11 (7,8)
Pipe#5 (3,6) Pipe#8 (5,6)
Loop #4
8
Pipe#12 (8,9)
Demand 2 MMSCFD Fig. 3. Network topology and loops of Case 1.
6
Demand 2 MMSCFD
Pipe#10 (6,9)
Pipe#6 (4,5)
7
Demand 2 MMSCFD
3
Loop #2
5 Pipe#9 (5,8)
4 Pipe#7 (4,7)
Demand 3 MMSCFD
Pipe#3 (2,3)
Pipe#4 (2,5)
1 Pipe#2 (1,4)
Supply 16 MMSCFD
9
Demand 2 MMSCFD psf = 130 psia
234
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Table 1 Summary of the loops in Case 1.
Linearized loop equations
Nodal mass balance equations
Fig. 4 shows the resulting system of equations generated by the linear-rate analog method in a matrix form. Twelve linear equations are solved to update twelve flow rates in the pipeline branches. The first eight equations are nodal mass balance equations while the other four equations are linearized loop equations. This system of equations can be solved utilizing any linear solver. For comparison purposes, three linear solvers are used to solve this problem: the linear-rate analog (this work), linear-pressure analog (Ayala and Leong, 2013), and the Linear theory method (Wood and Charles, 1972). For the linear theory solution, a random distributed gas flow rates between 1 and 2 MMSCF/day is used as initial flow guess because it is required for linear-theory iterations to begin. Table 2 and Table 3 show that the prediction of node pressures and flow rates are identical by implementing any of the three aforementioned methods. Fig. 5 shows that it takes 38 and 34 iterations and 6.705 ms and 6.399 ms of CPU time for the linear theory and linear-pressure analog, respectively, to achieve the prescribed flow rate tolerance (1 104 SCF/day). Linear-flow analog takes, in turn, only 19 iterations and 4.836 ms of CPU time to achieve the same convergence of 1 104 SCF/day. The linearrate analog also exhibits a more robust performance than linear theory for this case study. An important reason for this is that the resulting coefficient matrix is not explicitly dependent of flow rates in the linear-rate case, which significantly speeds up convergence and eliminates the need for flow-rate initializations. In addition, the linear-rate analog takes fewer iterations than even its linearpressure counterpart, for the same identical case study. This can be explained by comparing the non-linear complexity of both formulations: the system of equations generated by the q-formulation contains less nonlinear equations than that of p-formulation.
Table 2 Summary of node pressure prediction using different methods, Case Study 1. Nodal pressure, psia Node
Linear-rate analog
Linear-pressure analog (Ayala and Leong, 2013)
Linear theory
1 2 3 4 5 6 7 8 9
536.112 378.833 328.874 351.089 264.165 198.796 300.382 187.528 130
536.112 378.833 328.874 351.089 264.165 198.796 300.382 187.528 130
536.112 378.833 328.874 351.089 264.165 198.796 300.382 187.528 130
Table 3 Summary of flow rate predictions using different methods, Case Study 1. Gas rate, SCF/day Inlet
Outlet
Linear-rate analog
Linear-pressure analog (Ayala and Leong, 2013)
Linear theory
1 1 2 2 3 4 4 5 5 6 7 8
2 4 3 5 6 5 7 6 8 9 8 9
7 736 768.630 8 263 231.370 3 834 954.289 1 901 814.340 1 834 954.289 1 619 725.416 3 643 505.954 1 218 428.476 1 303 111.281 1 053 382.765 1 643 505.954 946 617.235
7 736 768.630 8 263 231.370 3 834 954.289 1 901 814.340 1 834 954.289 1 619 725.416 3 643 505.954 1 218 428.476 1 303 111.281 1 053 382.765 1 643 505.954 946 617.235
7 736 768.630 8 263 231.370 3 834 954.289 1 901 814.340 1 834 954.289 1 619 725.416 3 643 505.954 1 218 428.476 1 303 111.281 1 053 382.765 1 643 505.954 946 617.235
-1
-1
0
0
0
0
0
0
0
0
0
0
q1
-16000000
1
0
-1
-1
0
0
0
0
0
0
0
0
q2
2000000
0
0
1
0
-1
0
0
0
0
0
0
0
q3
2000000
0
1
0
0
0
-1
-1
0
0
0
0
0
q4
3000000
0
0
0
1
0
1
0
-1
-1
0
0
0
q5
1000000
0
0
0
0
1
0
0
1
0
-1
0
0
0
0
0
0
0
0
1
0
0
0
-1
0
q7
2000000
0
0
0
0
0
0
0
0
1
0
1
-1
q8
2000000
-λ1,4
0
λ2,5
0
- λ4,5
0
0
0
0
0
0
q9
0
0
0
λ2,3
- λ2,5
λ3,6
0
0
- λ5,6
0
0
0
0
q10
0
0
0
0
0
0
λ4,5
-λ4.7
0
λ5,8
0
- λ7,8
0
q11
0
0
0
0
0
0
0
0
λ5,6
- λ5,8
λ6,9
0
- λ8,9
q12
0
λ1,2
Fig. 4. Coefficient matrix of linear rate analog of Case Study 1.
×
q6
=
2000000
Fig. 5. Comparison of numerical performance of various linear solution techniques.
Fig. 6. Performance behavior - the analog linear-resistance transform t of the pipe branches.
236
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Fig. 7. Performance behavior - Gas flow rate solution.
Fig. 6, Fig. 7 and Fig. 8 display the convergence performance behavior of the linear-rate analog in terms of the evolution of the values of tij , flow rate, and pressure, respectively. As indicated in Fig. 2, the first iteration starts with tij ¼ to ¼ 1 for all pipes. It is noted in all these figures that convergence behavior of the proposed method is very smooth for all variables involved. Values of analog linear-resistance transform (tij ), flow rates, and pressure distributions become quite stable after the sixth iteration and updates for pipe branch flow rates and node pressures become small after few iterations. Additional numerical experiments were carried to quantify the impact of using different initial tij values (t0) on convergence behavior (other than the recommended tij ¼ 1 ¼ to used for the first iteration). As shown in Fig. 9, dramatically different initial t0 values (t0 < 0.01) were utilized compared to the to ¼ 1 choice. It is noted that convergence behavior of tij (and thus lij , since Cij is mostly or relatively constant depending on the choice of gas equation) is completely unaffected by the choice of to for the first iteration. This is due to the structure of the loop equations. For example, the expanded loop equation for Loop 1 (see Fig. 3) is shown below as Equation (17):
Equation (13), Equation (17) can be written as Equation (18):
t1;2
C1;2
q1
t1;4 C1;4
q2 þ
t2;5 C2;5
q4
t4;5 C4;5
q6
¼0
(18)
Therefore, for the very first iteration, this equation becomes:
t0
1 1 1 1 q q þ q q C1;2 1 C1;4 2 C2;5 4 C4;5 6
¼0
(19)
which demonstrates that using any value of t0 would not alter results from the first iteration given that such value cancels out of the equation (Equation (19)). Hence, unlike the linear-pressure analog where different initial linear-conductivity initializations can be shown to have some impact on convergence, it is demonstrated that convergence behavior of the proposed linear-rate analog is not sensitive to the t0 values used in the first iteration, and that the recommended use of lij ¼ 1=Cij (to ¼ 1) would suffice.
5.2. Case study 2: inclined pipe network system
l12 q1 l14 q2 þ l25 q4 l45 q6 ¼ 0
(17)
Expressing the linear rate analog resistivity terms using
In order to consider the elevation of the pipeline network, Equation (4) is modified written as Equation (20) (Ferguson, 1936):
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
237
Fig. 8. Performance behavior - Nodal pressure solution.
p2i esij p2j ¼
qij Cij
the inclined pipe Equation (20) can be rewritten for the case of a linear-q formulation as follows:
!2 (20)
where sij is a term which characterizes the elevation different between the upstream and downstream nodes, which is defined in Equation (21) (Ayala, 2013):
sij ¼ 0:0375
gg D H zave TAve
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 sij u " sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi # u pi e 2 pj 1 2 7 61 u q pi e pj ¼ 4 t ¼ 1 q 5 ij sij Cij C 1 þ rij ij ij 2 pi þ e pj 2
sij 2
¼
1 t Cij ij
(21)
(24)
In this case, the pipeline conductivity is written in terms of equivalent pipe lengths, as shown in Equation (22) (Ayala, 2013):
where the linear-rate analog resistance (lij Þ and the analog linearresistance transform (tij) are again defined for the inclined case as:
Ef sG pTscsc 1 d2:5 Cij ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi 0:5 gg Tave zave fF Le
1 lij ¼ Cij
(22)
In Equation (22), Le is the equivalent pipeline length, in ft, which can be calculated as (Ayala, 2013):
Le ¼
ðesij 1Þ L sij
(23)
Following the development presented in Equations (10)e(14),
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 1 ¼ t 1 þ rij Cij ij
(25)
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 tij ¼ 1 1 þ rij
(26) sij
but with rij ¼ spij i . Note the presence of the ‘e 2 ’ in denominator of 2 pj the rij ratio foreinclined pipes. Therefore, the linear-flow analog flow equation model proposed by this work for the inclined case is:
238
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Fig. 9. Impact of initial values transformation variable t on the iteration efficacy.
sij
pi e 2 pj ¼ lij qij
(27)
Or,
sij pi pj ¼ lij qij þ e 2 1 pj
(28)
As a result, Kirchhoff's second law (Equation (16)) is rewritten as follows for the inclined case as follows: loop loop X h sij i X ± pi pj ¼ xij lij qij þ e 2 1 pj ¼ 0 i;j
(29)
ij
or, loop X ij
xij lij qij ¼
loop X
sij
xij 1 e 2 pj
(30)
ij
Again, xij equals to positive one if the flow direction agrees with the prescribed loop positive direction, or equals to negative one if it goes against. The strategy is similar to the horizontal case, except that the value of the right-hand side (RHS) of Equation (30) needs to be updated at every iteration other than the first. For the first
iteration, tij ¼ 1 ¼ to and the RHS of Equation (30) is ignored (as if it were a horizontal system) to startup the iteration process. Linear-rate analog method is now implemented to solve the inclined pipeline network proposed by Ayala and Leong (2013), and displayed in Fig. 10. This network system is similar to Case study 1 (four loops and nine nodes) but nodes at located at different elevations. Results obtained via the implementation of the linear-rate analog scheme are presented in Table 4 and Table 5 in terms of nodal pressure and pipeline branch flow rates, respectively. It is noteworthy that the solution of linear rate analog shows an exact agreement with the reported results generated by Ayala and Leong using linear pressure analog. As shown in Fig. 11, however, the linear-rate analog takes 19 iterations to achieved the prescribed flow rate tolerance of 1 104 SCF/day, while linear-pressure analog takes 34 iterations and linear theory takes 39 iterations to achieve the same level of tolerance. Also, the required CPU time for linear-rate analog is smaller compared to the other linearization techniques. It should be added that the linear theory method was not able to originally converge for this case study until all known flow directions were explicitly specified (which was done in order to generate the comparison in Fig. 11). Neither linear-pressure nor linear-rate analogs require for any branch flow directions to be known to ensure convergence.
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
239
Demand 2 MMSCFD
7
Δh (7,8) = 600 ft
Pipe#11 (7,8)
Pipe#8 (5,6)
8
Loop #4
Δh (8,9) = -1000 ft
Pipe#12 (8,9)
Pipe#5 (3,6)
Δh (3,6) = - 800 ft
Δh (2,5) = - 400 ft
Pipe#4 (2,5)
Δh (5,6) = 1200 ft
6
Demand 2 MMSCFD
Pipe#10 (6,9)
Loop #3
Loop #2
Δh (6,9) = -1600 ft
Demand 1 MMSCFD
Demand 2 MMSCFD
3
Pipe#3 (2,3)
5
Pipe#6 (4,5) Pipe#7 (4,7)
Δh (4,7) = - 800 ft
Pipe#1 (1,2)
Loop #1
Δh (2,3) = 1200 ft
2
Δh (4,5) = -800 ft
4
Demand 2 MMSCFD
Δh (1,2) = -1200 ft
Pipe#9 (5,8)
Demand 3 MMSCFD
Pipe#2 (1,4)
Δh (1,4) = - 800 ft
1
Δh (5,8) = 600 ft
Supply 16 MMSCFD
9
Demand 2 MMSCFD psf = 130 psia
Demand 2 MMSCFD Fig. 10. Network topology and identified loops for inclined pipeline network system.
6. Concluding remarks
Table 4 Summary of the node pressure prediction, Case study 2. Node pressure, psia Node
Linear-rate analog
Linear-pressure analog (Ayala and Leong, 2013)
Linear theory
1 2 3 4 5 6 7 8 9
524.759 380.602 319.392 346.358 266.612 191.583 302.400 185.047 130
524.759 380.602 319.392 346.358 266.612 191.583 302.400 185.047 130
524.759 380.602 319.392 346.358 266.612 191.583 302.400 185.047 130
Table 5 Summary of the gas flow rate prediction, Case study 2. Gas rate, SCF/day Inlet
outlet
Linear-rate analog
Linear-pressure analog (Ayala and Leong, 2013)
Linear theory
1 1 2 2 3 4 4 5 5 6 7 8
2 4 3 5 6 5 7 6 8 9 8 9
7 8 3 1 1 1 3 1 1 1 1 9
7 8 3 1 1 1 3 1 1 1 1 9
7 8 3 1 1 1 3 1 1 1 1 9
742 346.142 257 653.858 810 225.731 932 120.412 810 225.731 612 367.344 645 286.513 231 468.278 313 019.478 041 694.009 645 286.513 58 305.991
742 346.142 257 653.858 810 225.731 932 120.412 810 225.731 612 367.344 645 286.513 231 468.278 313 019.478 041 694.009 645 286.513 58 305.991
742 346.142 257 653.858 810 225.731 932 120.412 810 225.731 612 367.344 645 286.513 231 468.278 313 019.478 041 694.009 645 286.513 58 305.991
In this work, a linear-rate analog formulation is proposed to solve the natural gas transportation in pipeline networks problems using a q- (flow-) formulation. Linear-rate analog effectively linearizes the loop equations formulated by Kirchhoff's second law of a gas pipeline network system via the use of linear-rate analog pipeline resistivities. Through case studies, the linear-rate analog is shown to be a robust, stable and initial guess-free solution technique for natural gas pipeline network analysis. The following additional observations may be drawn: Linear-rate analog can be successfully implemented to solve both horizontal and inclined pipeline network cases. The solutions of linear-rate analog agree well with known reported results in terms of node pressures and branch flow rates for both cases. The numerical performance of the linear-rate analog method has been compared against two other avalaible linearization algorithms, linear theory and linear-pressure analog method, by solving the same pipeline network cases. It is shown that the linear-rate analog requires fewer iteration for convergence and smaller CPU time than the other two solution techniques for all cases under consideration. The linear-rate analog method is an initial-guess free technique (no pressure or flow rate guesses are needed to start iterations) that linearizes the loop equations and constructs a coefficient matrix that is independent of primary unknowns (flow rate). Computational savings associated with the use of linear-rate analog method is expected to grow considerably as the size of the network system increases because most of the equations deployed by the q-formulation are linear equations.
240
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
Fig. 11. Comparison of numerical performances of linear-rate analog against linear theory and linear-pressure analog e Case Study 2.
Nomenclature AGA Cij Di d e fF DH L pi qij r rG R Rij s Si T z
American Gas Association pipeline conductivity of a pipe branch connected by nodes i and j SCF/day-psi gas demand at an arbitrary node i SCF/day pipeline diameter inch pipeline toughness ft Fanning friction factor dimensionless elevation difference between nodes ft pipeline length ft pressure of an arbitrary node i psia gas flow rate of a pipe branch connected by nodes i and j SCF/day pipeline node pressure ratio dimensionless pipeline specific resistivity psia2-inch5/SCFD2-ft gas constant, with a magnitude of 10.732 psia-ft3lbmole1- R1 pipeline resistivity of pipe branch connected by node i and j psia2/SCF2/day2 pipe elevation parameter dimensionless gas supply at an arbitrary node i SCF/day temperature R gas compressibility factor dimensionless
Appendix A. Summary of gas flow equations
Table A-1 Summary of the gas equations (Modified from Ayala, 2013). Generalized gas equation in conduits: qij ¼ Cij Ef sG
Tsc psc
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2i esij p2j ,
Cij ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p1ffiffiffi gg Tave zave
d2:5 0:5 fF Le
Gas equation General equation
Definition of friction factor ! p1ffiffiffi ¼ 4:0 log10 fF
Weymouth
fF ¼ 0:008 d1=3
Panhandle-A
fF ¼
0:01923
fF ¼
0:00359
qgSC gg d
Panhandle-B
qgSC gg d
AGA Fully Turbulent
3:7 þ
0:03922
3:7d e
where Re is the Reynolds number, qgSC is in SCF/day, L is in ft or mile, d is in inch, p is in psia and T is in R. The unit conversion factor sG varies based on the unit of the pipe length. For L in ft, sG ¼ 2818, sG,w ¼ 31 508, sG,PA ¼ 44 400, sG,PA ¼ 58 328; For L in miles, sG ¼ 3 ¼ G,w ¼ 4 ¼ sG,PA ¼ 435.98, sG,PA ¼ 736.
Greek
g l s x t
gas specific gravity dimensionless Linear-rate analog pipeline resistivity psi/SCF/day gas equation unit conversion factor directional sign term of a close pipeline loop transformation variable
Subscripts ave average value e equivalent i branch upstream node j branch downstream node sc standard condition 0 first iteration
5:02 pffiffiffi Re fF
0:1461
p1ffiffiffi ¼ 4 log10 fF
e d
Appendix B. Network details e Case Studies
Table B-1 Gas properties and general data (Case Study 1 and 2). Base pressure (SC)
14.7 psia
Base temperature (SC) Average gas compressibility Average gas gravity T average Pipe roughness Pipe efficiency
60.0 F 0.9 0.69 535 R 0.0018 in 1
Q. Sun, L.F. Ayala / Journal of Natural Gas Science and Engineering 43 (2017) 230e241
241
Table B-2 Pipe and supply/demand data (Case Study 1 and 2).
References American Gas Association AGA, 1965. Steady Flow in Gas Pipelines: Testing,/ Measuremen,/Behavior, and Computation. Institute of Gas Technology, New York, NY, U.S. Ayala, H.,L.F., 2013. Chapter 22: Transportation of crude oil, natural gas, and petroleum products. In: Petroleum Refining and Natural Gas Processing. ASTM International, West Conshohocken, PA, USA, pp. 549e575. Ayala, L.F., Leong, C.Y., 2013. A robust linear-pressure analog for the analysis of a natural gas transportation networks. J. Nat. Gas Sci. Eng. 174e184. Brki c, D., 2009. An improvement of Hardy Cross method applied on looped spatial. Appl. Energy 86, 1290e1300. Chen, N.H., 1979. An explicit equation for friction factor in pipe. Ind. Eng. Chem. Fund. 18, 296e297. Colebrook, C.F., 1939. Turbulent flow in pipes, with particular reference to the transition region between smooth and rough pipe laws. J. Inst. Civ. Eng. 12, 133e156. Cross, H., 1930. Analysis of continuous frames by distributing fixed-end moments. Proc. Am. Soc. Civ. Eng. 56 (5), 919e928. Energy Information Administration EIA, 2013. Technically Recoverable Shale Oil and Shale Gas Resources: an Assessment of 137 Shale Formations in 41 Countries
outside the United States. EIA Independent Statistics and Analysis. U.S. Energy Information Administration, Washington DC. U.S. Ferguson, J.W., 1936. A mathematical analysis of the effect of differences on flow formulae for gas pipe lines. Gas. Age 78, 72e76. Menon, E.S., 2005. Gas Pipeline Hydraulics. Taylor and Francis Group, Boca Raton, FL, U.S. Gas Processors Supplier Association GPSA, 2004. “Fluid Flow and Piping.” GPSA Engineering Data Book. 12. Gas Processors Supplier Association, Tulsa, OK, U.S. Kelkar, M., 2008. Natural Gas Production Engineering. Penn Well Corporation, Tulsa, OK, U.S. Kumar, S., 1987. Gas production Engineering. Gulf Pub. Co., Houston. Leong, C.Y., Ayala, L.F., 2014. Hybrid approach using linear analogs for gas network simulation with multiple components. Oil Gas Facil. J. 3, 76e88. Makogon, Y.F., Holditch, S.A., Makogon, T.Y., 2007. Natural gas-hydrates-a potential energy source. J. Pet. Sci. Eng. 56, 14e31. Moody, L.F., 1944. Friction factors for pipe flow. Trans. ASME 66, 671e684. Weymouth, T.R., 1912. Problems in natural gas engineering. Trans. ASME 34, 185e206. Wood, D.A., Charles, C.O., 1972. Hydraulic network analysis using linear theory. ASCE J. Hydraul. Div. 1157e1170. Wood, D.A., Mokhatab, S., Economides, M.J., 2008. Technology options for securing markets for remote gas. In: Proceedings of 87th Annual Convention of the Gas Processors Association. GPA, Grapevine, Texas, US.