Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 41, 1, pp. 33-54, Warsaw 2003 FLOW MODELLING IN GAS TRANSMISSION NETWORKS. Part II – METHOD OF SOLUTION AND ALGORITHMS Robert Fedorowicz Edward Kołodziński Lech Solarz Military University of Technology, Warsaw e-mail: rfedor@ias.wat.waw.pl; solarz@ias.wat.waw.pl In the first part, the mathematical problem of flow modelling in gas transmission network was formulated. In the paper, the method of so- lution is explained for the model with the postulate T = const. The method of automatic generation of the system of algebraic equations is described. The results, good and bad experiences are described. Key words: flow in nets, gas transmissionmodelling, gas flow 1. Solution method The system of equations, boundary and initial conditions in dimensionless form are described for the isothermal model. 1.1. The problem in dimensionless form All dimensional quantities are related to characteristic quantities. For example: pref = 5MPa, Vref = 5m/s, tref = 3600s xref = 1000m, href = 1m, Tref equal to average temperature of the gas, ρref is calcu- lated on the basis of formula (3.6)4 in Fedorowicz et al. (2002) (Part I). pref = ρrefZ(ρref,Tref), qref = ρrefVref ζ = x xref τ = t tref (1.1) As the unknowns we take ρ(x,t) and dimensionless density of gas flux q(x,t) q= ρV = Q AρrefVref (1.2) 34 R. Fedorowicz et al. Wehave also used other pairs of the unknowns, but the pair ρ, q seems to be the most suitable. The equations of the problem are derived from (3.6)1,2 (Fedorowicz et al., 2002, Part I) ∂ρ ∂τ + 1 f1 ∂q ∂ζ =0 (1.3) ∂q ∂τ + 1 f1 ∂ ∂ζ (Vq)+ 1 f1f2 ∂ ∂ζ ( ρZ(ρ,1) ) + r3 f1 Vq+ρ r4 f1 dh(ζ) dζ =0 where f1 = xref trefVref f2 = (Vref aref )2 aref = √ pref[MPa] ·10 6 ρref r3 = xrefλ(Re,ε) 2D Re= q µN µ(ρ,1) Dr5 r5 = ρrefVref µN r4 = ghref V 2ref Z(ρ,1)= Z(ρρref,Tref) Z(ρN,TN) (1.4) The dimensionless pressure and temperature are equal to p= ρZ(ρ,1) T =1 (1.5) 1.2. Reduction of the partial differential equations to the algebraic equations Each pipe of the network is divided into Nip equal parts, Nip ∈ [1,64], ip is the index of the pipe ∆ζip = Lip xrefNip (1.6) The number Nip is chosen in such way that ∆ζip is almost the same for each pipe ip if it is possible. The time step is ∆τ = 1/Nτ , where Nτ is a natural number. The values of the unknowns at point ζk = k∆ζip at time τ = t∆τ, are ρk,t, qk,t, where k ∈ [0,1,2, ...,Nip], t ∈ [0,1,2, ...,Nτ ]. The unknown functions ρ, q, and V are approximated by straight lines between the points ζk. Integrating equations (1.3) in the intervals (ζk,ζk+1), we get Flow modelling in gas transmission networks... 35 dρk dτ + dρk+1 dτ +qk+1 2 f1∆ζip −qk 2 f1∆ζip =0 (1.7) dqk dτ + dqk+1 dτ + 2 f1∆ζip (Vk+1qk+1−Vkqk)+ + 2 f1∆ζipf2 [ρk+1Z(ρk+1,1)−ρkZ(ρk,1)]+ r3 ( Reav,ρav) 3f1 · ·[(2Vk +Vk+1)qk+(Vk+2Vk+1)qk+1]+(ρk+ρk+1) r4 f1 dh(ζ) dζ =0 The symbol av is explained under the formula (1.11). The τ derivatives are approximated by the difference quotients dρk dτ = 1 ∆τ (ρk,t−ρk,t−1) dρk+1 dτ = 1 ∆τ (ρk+1,t−ρk+1,t−1) (1.8) but in the remaining terms of equations (1.7), the unknowns ρk, qk are repla- ced by ρk ⇒αρk,t+(1−α)ρk,t−1 qk ⇒αqk,t+(1−α)qk,t−1 (1.9) Unknowns ρk+1, qk+1 are replaced in the sameway. Theweight α∈ (0,1). The influence of the value of α, on the results of calculations, was investigated in Fedorowicz (2001). When α=0.5, the method is equivalent to the Crank- Nicholson method for linear equations. Rearranged equation (1.7)1 assumes the form qk,t(−wip)+ρk,t+qk+1,twip+ρk+1,t = (1.10) = ρk,t−1+ρk+1,t−1+νip(qk,t−1− qk+1,t−1) where wip = 2∆τ f1∆ζip α νip = 2∆τ f1∆ζip (1−α) Rearranged equation (1.7)2 takes the form qk,t(1+wipBk)+ρk,twipAk+qk+1,t(1+wipBk+1)+ρk+1,twipAk+1 = =−νip(ρk,t−1Ak+qk,t−1Bk+ρk+1,t−1Ak+1+qk+1,t−1Bk+1)+ (1.11) +qk,t−1+qk+1,t−1 where Ak = r4∆ζip 2 dh dζ − 1 f2 Z(αρk,t+(1−α)ρk,t−1,1) 36 R. Fedorowicz et al. Ak+1 = r4∆ζip 2 dh dζ + 1 f2 Z(αρk+1,t+(1−α)ρk+1,t−1,1) Bk =−Vk+ ∆ζip 3 r3(Reav,ρav) ( Vk+ 1 2 Vk+1 ) Bk+1 =Vk+1+ ∆ζip 3 r3(Reav,ρav) ( Vk+1+ 1 2 Vk ) Vk = αqk,t+(1−α)qk,t−1 αρk,t+(1−α)ρk,t−1 Vk+1 = αqk+1,t+(1−α)qk+1,t−1 αρk+1,t+(1−α)ρk+1,t−1 qav = 1 2 {[αqk,t+(1−α)qk,t−1]+ [αqk+1,t+(1−α)qk+1,t−1]} ρav = 1 2 {[αρk,t+(1−α)ρk,t−1]+ [αρk+1,t+(1−α)ρk+1,t−1]} The values qav, ρav are used for calculation of Reav. 1.3. Reordering of the unknowns Theunknowns ρk,t, qk,t are ordered in the vector x(No), where the symbol No denotes the consecutive number of the unknown. The rule of ordering ρk, qk =⇒x(m) is illustrated in Fig.1 Mip =2 [ ip ∑ i=1 (Ni+1)− (Nip+1) ] (1.12) Fig. 1. The rule of ordering Flow modelling in gas transmission networks... 37 Taking into account the equations (1.10), (1.11), boundary conditions de- scribed inPart I, and the rule of ordering ρk, qk =⇒x(m) the software packet generates the system of non-linear algebraic equations J ∑ j=1 Aij(x1,x2, ...,xJ,xp1,xp2, ...xpJ)xj =Ri(x1,x2, ...,xJ,xp1,xp2, ...xpJ, t) (1.13) i=1,2, ...,J J =2 IP ∑ ip=1 (Nip+1) where IP is the total number of pipes in the connected subnetwork.Whenwe calculate {−→x}, the values of ρk,t−1 and qk,t−1 are known and ordered in the vector −→xp. 1.4. Iterative linearization of the nonlinear algebraic equations The system of nonlinear algebraic equations (1.13) is reduced to the con- secutive systems of linear algebraic equations J ∑ j=1 Aij ( x (k−1) 1 ,x (k−1) 2 , ...,x (k−1) J ,xp1,xp2, ...,xpJ ) x (k) j = (1.14) =Ri ( x (k−1) 1 ,x (k−1) 2 , ...,x (k−1) J ,xp1,xp2, ...,xpJ, t ) where i=1,2, ...,J, k=1,2, ...,K. For k = 0 we take x (0) j = xpj, j = 1,2, ...,J. The matrix A = {Aij} is sparse. Themeasures ∆ρσ = √ √ √ √ √ √ J/2 ∑ i=1 (x (k) 2i −x (k−1) 2i ) 2 J/2 ∆qσ = √ √ √ √ √ √ J/2−1 ∑ i=0 (x (k) 2i+1−x (k−1) 2i+1 ) 2 J/2 (1.15) are used to investigate the convergence.The measures (1.15) are related to ρσ, qσ ρσ = √ √ √ √ √ √ J/2 ∑ i=1 (x (0) 2i ) 2 J/2 qσ = √ √ √ √ √ √ J/2−1 ∑ i=0 (x (0) 2i+1) 2 J/2 (1.16) The fact that the iterationprocess is divergent, is signalizedby theprogram and calculations in this subnetwork are stopped. This happenswhenwe try to 38 R. Fedorowicz et al. simulate the transportof a too large amountof gas, butalso,whenwecalculate the flow in the middle of the summer, so, when the flow ratio is very small. We can influence the speed of convergence by the changes of ∆τ and ∆ζip. The changes of ∆ρσ/ρσ and ∆qσ/qσ as functions of kwere investigated. The qmeasure ismore susceptible to the fluctuation than the ρmeasure. After six or seven iterations the measures of distance between iterations remain rather constant. Our customers prefer the calculations with very small numbers of iterations (K =1or K =2) because the precision of calculation is practically not very important for them. The customers knowledge of the parameters of the network is rather small, so the precision of the mathematical solution is very high from their point of view. 1.5. Solution of the system of linearized algebraic equations One of the direct methods is applied – the Gauss elimination with back substitution and row pivoting. The software realization is different than in Pissanetzky (1984), Duff et al. (1986), because it is based on ordered storage of the A matrix in the list. It seems to be a better solution when the class ”list” exists. The single linked list is used. The interesting details are given in the Appendix. 2. Generation of the nonlinear system of algebraic equations 2.1. Gaining the information about the network, the gas, the network control and load 2.1.1. Gaining the information about the network The description of practical aspects of gaining information about real transmission networks is contained in the conclusions. This information is in principle imported from the database of GIS type. The packet has the interfa- ce which enables the conversion and edition of the data. The usermay change the data, and even hemay introduce a completely new network. Each pipe is characterized by diameter D, length L, relative roughness ε, number N. Each pipe is identified by its registration numbers. The first and the second nodes incident with the pipe are identified by their registration number. The consecutive number ip and the Boolean variable ”hasbeen” are also included in the CPipe class. Each node is characterized by the height h, the natural number ”whatnode” for the identification of the node type, the coefficient of local Flow modelling in gas transmission networks... 39 pressure loss Λ, and the registration number of the incident directional pipe (see Fig.6, Part I). The node is identified by its registration number. In the CNode class we findmemberswhich characterize the one-sided constrains (see Subsection 2.1, Part I). As in the CPipe class, the consecutive number and the Boolean variable ”hasbeen” are also included in the CNode class. The information about the network is stored in the tables PIPES, NODES. 2.1.2. Gaining the information about the gas The data about the gas emerge from the chromatographic analysis. The packet enables including the data into the small database, which is a part of the packet. User-friendly interface enables the change of compounds, mixing the new gas with the gases existing in the base. For the flow calculation, only some parameters which are enumerated in Subsection 2.2 of the Part I are necessary. Additionally, we include the variables which characterize the type of the gas (high methane, high nitrogen). The type of the gas influences the formulae µN/µ(ρ,1) and Z(ρ,1), (1.4). 2.1.3. Information about control of the network The external nodes, such as pressure sources, sources with known flow-in ratios, whichmay be controlled, get the control functions. The internal nodes with compressor stations, pressure reducers, valves, gate valves also get the control functions. The control function is a continuous step function, such as that shown in Fig.2. Fig. 2. Control function 40 R. Fedorowicz et al. The data concerning the control function are stored in the table OPERATION∈CCF class. The members of CCF class are: the type of the node ”whatnode”, the registration number of the node, the mode of control (pressure, flow ratio, compression), the number of pairs of the data (ti,CFi), and the table with the pairs. The table OPERATION is common for all sub- networks. 2.1.4. Information about the loads The information about the consumption of gas through the reduction sta- tions of the first step (city gate), about certain sources of gas, are forecasted. These data are not controlled. Each element of these data contains the daily load Qd ormean pressure for the pressure source, and the consecutive number of the profile in the base of different profiles (see 2.4, Part I). Taking these data into account, we prepare the table of data, named SFORECAST, about dimensionless forecasted: flow-off ratio, flow-in ratio, pressure. Each element of the table contains: the registration number of the node, and tables ai, bi, ci, Yi; i=1,2, ...,48. The tables determine coefficients of the spline interpolation of the data Y (τ)= [(aiξ+bi)ξ+ ci]ξ+Yi ξ= τ − int(τ) (2.1) where the function int(·) returns the integer part of the argument. The function (2.1) is determined for τ ∈ [0,48] but it is used only for 24 hours. The table SFORECAST is common for all subnetworks. 2.2. Gaining the initial conditions The flow calculations simulating the process remaining constant in time (static model) are enabled in the packet. The static calculations are the only possibility of gaining information about the flow in the initial moment of cal- culation, in certain cases. The method of static calculations is not described in this paper. The last results of the previous calculations can be taken as the initial conditions. We shortly name such situation the continued calculations. The data are of the same standard,with the registration numbers of pipes included into the data. It is a flexible system. For example, the results of calculations of a more robust network may be taken as the initial conditions for a smaller network. The initial conditions are stored in −→xp. Flow modelling in gas transmission networks... 41 2.3. Separation of connected subnetworks 2.3.1. Closed valves The data concerning the network are given in tables NODES, PIPES.We search the closed valves and gate valves in NODES and we check the table OPERATION. If the valves are closed, then the virtual nodes are added with the same registration numbers the as the existing nodes, but with the sign ”minus” before them. For all these nodes the flow ratio =0. 2.3.2. Information about a subnetwork. The table ”incy” When building equations (1.13), we apply the information about the con- nected subnetwork (the subnetwork with connected graph) stored in a two- dimensional table ”incy”. The name of the table suggests likeness to an inci- dence matrix, but it contains much more data. The table enables us to find all the necessary data once, not at all steps t. The second index of the table is for the consecutive number of the pipe in ordered table. The meaning of the first index is illustrated in Fig.3 Fig. 3. The table ”incy” Below, the symbols =⇒ and ⇐= are used. The symbol =⇒ means that the considered pipe has the same direction as the pipe ip, the symbol ⇐= is for the opposite direction. The consecutive position of the table are for: • incy(1, ip), the consecutive number. In the final stage equal to ip • incy(2, ip), the consecutive number of the first pipe incident with the second node of the ip-pipe; if there is not such a pipe, then incy(2, ip) = 0; if =⇒ then incy(2, ip)> 0; if ⇐=then incy(2, ip)< 0 • incy(3, ip), the consecutive number of the second pipe incident with the second node of the ip-pipe; if there is no second pipe then 42 R. Fedorowicz et al. incy(3, ip) = 0; if incy(2, ip) = 0 then incy(3, ip) = 0; always |incy(3, ip)|> |incy(2, ip)| if incy(3, ip) 6=0, if =⇒ then incy(3, ip)> 0; if ⇐= then incy(3, ip)< 0 • incy(4, ip), the consecutive number of the first pipe incident with the first node of the ip-pipe; if there is no such pipe then incy(4, ip) = 0; if =⇒ then incy(4, ip)> 0; if ⇐= then incy(4, ip)< 0 • incy(5, ip), the registration number of the first node of the ip-pipe • incy(6, ip), the registration number of the second node of the ip-pipe • incy(7, ip), the consecutive number in the table SFORECAST conta- ining data about the loads or flow-in ratio or the consecutive number (plus the bias number) in the table OPERATION. The number is po- sitive. If incy(7, ip) = 0 then the first node of the ip-pipe is without known loads or flow-in ratio. If the node is a compression station and the ip-pipe is not its directional pipe, then also incy(7, ip) = 0. If the incy(7, ip) contains the consecutive number in the table SFORECAST, then it is less than in the case of the table OPERATION • incy(8, ip), the consecutive number of the first node of the ip-pipe in the ordered table of nodes • incy(9, ip), the consecutive number of the second node of the ip-pipe in the ordered table of nodes • incy(10, ip), the consecutive numberof the secondpipe incidentwith the firstnode of the ip-pipe; if there is no secondpipe, then incy(10, ip) = 0; if incy(4, ip) = 0 then incy(10, ip) = 0; always |incy(10, ip)| > |incy(4, ip)| if incy(4, ip) 6= 0, if =⇒ then incy(10, ip) > 0; if ⇐= then incy(10, ip)< 0 • incy(11, ip) contains the consecutive number meaning the same as incy(7, ip) but for the second node of the ip-pipe. 2.3.3. Separation of the connected subnetworks from the network ”Combing”.Wecomb the tablesNODES andPIPES and take all nodes and pipes with the consecutive number equal to 0 to the tables NODESNEWand PIPESNEW.Thesenodesandpipeshavenotbeen included into thepreviously ”combed” subnetworks. If we ”comb” less than two nodes or no pipe, then the process of separation is finished. Flow modelling in gas transmission networks... 43 The node number one. We select the node with the consecutive number equal to 1 from the nodes in the table NODESNEW. The criteria: • The nodes – pressure sources or the nodes – gasholder with compressor stations in ”pressure mode”. From all such nodes the node with the highest pressure is chosen. If such nodes do not exist, then • The nodes – sources with known flow-in ratio or the nodes – gasholders with compressor stations in ”flow-off mode”. From all such nodes, the node with the highest absolute value of the flow ratio is chosen. If such nodes do not exist, then • The nodes – closed valves or gate valves. From all such nodes the node incident with the pipewith the highest diameter is chosen. If such nodes do not exist, then • The nodes – external nodes with forecasted flow-off ratio. From all such nodes, thenodewith thehighestabsolutevalueof theflowratio is chosen. If such nodes do not exist, then • The flow in the subnetwork without any loads or supplies of gas is not calculated. How, without any valves or gate valves, the gas has been brought in? The transient calculations of flow are enabled in the subnetwork without the source of gas of a known pressure. It is an important difference between static and transient models. The subnetwork ordering. The node number 1 is stored in the table NO- DESORDERED as number 1. The incident pipe gets number 1 in the table PIPESORDERED.The fact that a pipe or a node is stored in the ordered ta- bles is given in the tables NODES andPIPESwhere the consecutive numbers are changed from 0. The service of the first node and the first pipe is different from the service of the remaining pipes and nodes. If the node no.1 is the first node of the pipe no.1 then the positions 1,4,10,5,8,7,6,9, of incy(∗,1) get the values. If the node no.1 is the second node of the pipe no.1 then the positions 1,2,3,11,5,8,7,6,9, of incy(∗,1) get the values. The data about the adjacent node are stored in NODESORDERD[2]. Themembers ”hasbeen” of the node no.1 (and all consecutive nodes stored in NODESORDERED) get the value TRUE. The service of the node no.2.We count and register the incident pipes. 44 R. Fedorowicz et al. Case 1. If there is 1 incident pipe then the position 2,3,11 of incy(∗,1) get the values, when the node no.1 is the first or the positions 4,7,10 of incy(∗,1) get the values when the node no.1 is the second node of the pipe no.1.We define the structure as a tree andwe stop the subnetwork analysis, because it consist of one pipe and two nodes. Case 2. If there are 2 incident pipes then the positions 2,3,11 or 4,7,10 of incy(∗,1) get the values as in the previous case, so the service of the pipe no.1 is finished. The data about new incident pipe are stored in PIPESORDERED[2]. The data about new adjacent node are stored in NODESORDERD[3]. We begin the service of the incy(∗,2). The posi- tions 1,4,10,5,6,8,9, also 7, get thevalues if thenodeno.2 is thefirstnode of the pipe no.2 or the positions 1,2,3,5,6,8,9, also 11, get the values if the node no.2 is the second node of the pipe no.2. When we bind the nodewith the elements of the table SFORECASTorOPERATION, the ”directional pipe” is taken into account. Case 3. If there are 3 incident pipes then the diameters of the new pipes are investigated.Thenewpipewith grater diameter gets thenumber2and is stored inPIPESORDERED[2].Thenodeadjacent through this pipegets no.3 and is stored in NODESORDERD[3]. The new pipe with smaller diameter gets the number 3 and is stored in PIPESORDERED[3]. The node adjacent through this pipe gets no.4 and is stored in NODESOR- DERED[4]. The node no.2 is a three-way node. The positions 2,3,11 or 4,7,10 of incy(∗,1) get the values as in the first case, so the service of the pipe no.1 is finished. We begin the service of pipes no.2 and no.3. The positions 1,4,10,5,6,8,9, also 7, or the positions 1,2,3,5,6,8,9, also 11, get the values. After this analysis we have ”Jnod” nodes in the table NODESORDERED and”Ipip”pipes in the tablePIPESORDERED.Webegin the service of nodes from no.3 and we finish when, after the analysis, the consecutive number of the node is equal to ”Jnod”. At each step of the node service the nodes and pipes may be stored in the ordered tables and the value of ”Jnod” may grow up. For each analyzed node, we count and register the incident pipes. Case 4. If there is 1 incident pipe, the positions 2,3,11 or 4,7,10 of this inci- dent pipe get the values. Case 5. If there are 2 incident pipes, then we investigate if there is one of themwhich has not appeared in the table PIPESORDERED. Flow modelling in gas transmission networks... 45 1. If there is a new pipe then: we add 1 to ”Ipip”, we store the data of the pipe in PIPESORDERED[Ipip], if the node adjacent thro- ugh this pipe has not existed in NODESORDERED then we add 1 to ”Jnod” and we store the data about this node in the table NODESORDERED[Jnod]. The positions 2,3,11 or 4,7,10 of incy, for the second of the incident pipe, get the values. The positions 1,4,10,5,6,8,9, also 7 or 1,2,3,5,6,8,9 also 11 of the incy(∗,Ipip) get the values. 2. If there is no new pipe, then the positions 2,3,11 or 4,7,10 of incy, for the both incident pipes, get the values. Case 6. If there are 3 incident pipes then we investigate which of them have existed in the table PIPESORDERED. 1. If all of the abovementioned pipes have existed in the table PIPE- SORDERED, then the positions 2,3,11 or 4,7,10 of incy, for these incident pipes, get the values. 2. If there is one new pipe then: we add 1 to ”Ipip”, we store the data of the pipe in PIPESORDERED[Ipip], if the node adjacent through this pipe has not existed in NODESORDERED, then we add 1 to ”Jnod” and we store the data about this node in the table NODESORDERED[Jnod]. The positions 2,3,11 or 4,7,10 of incy, for the remaining incident pipes, get the values. The positions 1,4,10,5,6,8,9, also 7 or 1,2,3,5,6,8,9 also 11 of the incy(∗,Ipip) get the values. 3. If there are 2 new incident pipes then the diameters of the new pipes are investigated. We add 1 to ”Ipip”. The new pi- pe with greater diameter gets the number Ipip and is stored in PIPESORDERED[Ipip]. If the node adjacent to the investiga- ted node through this pipe has not existed in the table NODE- SORDERED, then we add 1 to Jnod and we store this node in NODESORDERED[Jnod].We add 1 to ”Ipip”. The newpipewith smaller diameter is stored in PIPESORDERED[Ipip]. If the no- de adjacent through this pipe did not exist in the table NODE- SORDERED, then we add 1 to Jnod and we store this node in NODESORDERED[Jnod]. The positions 2,3,11 or 4,7,10 of incy bound to the pipewhich exists inPIPESORDEREDget the values. We begin the service of newpipes.Thepositions 1,4,10,5,6,8,9, also 7 or the positions 1,2,3,5,6,8,9, also 11, get the values. 46 R. Fedorowicz et al. When this process is finished, then in Ipip elements of the tables incy, PIPESORDEREDand in Jnod elements of the table NODESORDERED,we have ordered data about the connected subnetwork. These subnetworks are counted and stored. 2.4. Generation of the system of algebraic equations 2.4.1. General remarks Wedonotpresent these excerpts of thealgorithmwhichconcern the testing and signaling incorrectness. An example of the incorrectness – the flow in the forbidden direction through the compressor station. The calculations are performed at the time t in all subnetworks. The cal- culations increase in time up to the moment when in all subnetworks the calculations are impossible, or up to the prescribed time. The subnetworks are calculated one after another in the same order as they were called into being (combed), so the calculations begin in the subnetwork with the highest pres- sure or the highest flow ratio. When the calculation for a time t are finished then the results are stored in the form suitable for visualizing themodule, and values from the table {−→x} are sent to the table {−→xp}. The storage of the sparse matrix {Aij} is described in the Appendix. 2.4.2. Generation of the system of algebraic equations for the subnetwork The generation is rather simple. All nodes in the NODESORDERED get the attribute hasbeen = FALSE. In the loop ip = 1,2, ...,IP, the actions described below are performed. • If the first node of the ip-pipe has the attribute hasbeen equal to FALSE, then thebelow enumeratedactions areperformed.Theattribu- te whatnode is analyzed. We include somany equations emerging from the boundary conditions as it is proper for the whatnode. The equations are coded in suchway, that the coefficient A of the xwith lower index is before the coefficient of the xwith greater index. The equations concer- ning the mass conservation principle are divided by the cross-sectional area of the ip-pipe, so we have the coefficient equal to 1. For the nodes with compressor station, reducer or check valve, the directional pipe is checked. We store the results in the list A, in the table of pointers p which point to the first coefficient in the equation (in the row), and in the table R containing the values of right side of the equation (see the Appendix). We count how many equations have been included and the first node of the ip-pipe gets the attribute hasbeen=TRUE. Flow modelling in gas transmission networks... 47 • We include Nip times equations emerging from (1.10) and (1.11) begin- ning from the nearest point to the first node of the ip-pipe. The rules concerning growing indices are applied. • If the second node of the ip-pipe has the attribute hasbeen equal to FALSE, then the analogues actions are performed as for the first node. As the result, we get the sparsematrix Aij with the elements concentrated near the main diagonal with splashes caused by boundary conditions. The preference of the pipes with greater diameters causes that the equations with more important unknowns are in the upper part of the matrix. 3. On the solution of the problem with heat exchange We do not take into account the dissipation of heat along the pipes. The explanation is in 4.1 (Part I) Fedorowicz et al. (2002). Two versions of the problem are possible: • thewhole network is calculatedwith the principle of energy conservation taken into account • the pipes near the compressor stations are calculated with the principle of energy conservation, but other pipes are calculated with the postulate T known. Submarine pipelines, pipelines in arctic zones (permafrost, marsh) must bemodelledwith the principle of energy conservation, but calculation ofmost pipelines in Polandmay be performed as it is described in the second section. We propose the method of solution shown below. • All dimensional quantities are related to characteristic quantities as in 1.1. • The partial differential equations are reduced to the algebraic equations as in 1.2. We get the system of equations d dτ ∑ i wlizi = ∑ m flm(z1,z2, ...,zM)zm (3.1) where l=1,2, ...,M, M – number of unknowns. The sumwith index i has two elements, the sum with index m has six, or four elements or even less. 48 R. Fedorowicz et al. • As in 1.2, the τ derivatives are approximated by difference quotients (1.8) and the unknowns ρk, qk,Tk are replaced by the sumswithweight α as in (1.9). We get the nonlinear set of algebraic equations analogous to (1.10), (1.11) but with the equation of energy included. • We apply the physical iteration. 1. The system of equations which emerged from the mass and mo- mentum equations, is solved with the unknowns Tk taken from the previous iteration. 2. The equations, which emerged from the energy equation, are solved with the above calculated unknowns ρk, qk. 3. The process described in points 1 and 2 is repeated and the conver- gence is tested. We start the iteration assuming the exponential decrease of temperature from the value known at the outlet of compressor station. The physical itera- tion process was tested for a very simple pipeline system with a compressor station (Komarec, 1999). The convergence was measured as in (1.15), (1.16). The convergence is rather fast, nearly one decimal point after each iteration (see also [8], Part I). 4. Results and predicted research The packet is ready and in use ([1], Part I). The calculations concerning different aspects of safety were performed, for example the survival time of a subnetworkwithout supply of gas, influence of the acceptance of new loads at the pressure in the subnetwork ([1, 5], Part I). The whole process of research and development of the packet yields the information about: simulation of the flow, methods of solution and the relations between the parameters we know and the parameters we should know. 4.1. Application of the modelling packet The accomplished research should prove to be useful. The packet is useful if we knowproper values of parameters and functionswhich influence the flow. Costs of usefulmeasurements are rather great, so this costs form the principal barrier whenwe think about broad application of anymodelling packet. If the packet is linkedwith thecomputerized systemof information then its recordsof Flow modelling in gas transmission networks... 49 repairs, inspections enable obtaining the necessary information ratherwithout costs, andvice versa, the computations enable verification of the data included into the system. The most difficult to obtain is the roughness of the pipes. The relations between: thematerial of the pipe, production technology, the time the pipe is in use, enable rough estimation of this parameter. The necessary next step is tuning of the roughness. Themeasured pressure and flow ratio are compared with the calculated values of pressure and flow ratio for different roughness of the pipes, and the best fitting roughness is taken as the proper value. Such process is realized even for the pipelines with very well known measured ro- ughness, e.g. the pipeline from Sleipner field to Zeebrugge. As a positive fact, we note that the roughness emerging from tuning in winter conditions fit the flow calculation in summer ([5], Part I). 4.2. Remarks about the methods of solution We get two successful methods of solution after some unsuccessful trials. 4.2.1. Uselessness of explicit finite differences schemes It was known that the explicit finite differences schemes are useless for this modelling, but the attempt to use this scheme to the 0 step of iteration seemed reasonable. When the iteration process begins, we take x (0) j = xpj, j=1,2, ...,J, so we take ∂ρ/∂τ =0 and ∂q/∂τ =0. It seems that if we take the time derivatives from the previous time step and use them to calculate x (0) j , we get better starting values for the iteration. The result obtained may be useful only when we show to students of the numerical methods, that the explicit methods are very doubtful. The figure of the function q(ζ,τ = const) looks like a saw blade. 4.2.2. Other useless methods The finite difference scheme of greater precision (central scheme) gave the system of equations for which the iteration process was very often divergent. Additionally, this central scheme of finite differences applied to the first-order differential equations leads to the second-order system of finite differences, so we had to add the new boundary conditions. We tried to solve the problem by another version of physical iteration.We have calculated the mass density ρ from the equations of mass conservation principle andnext, the flow ratio q from themomentumequation.Theprocess was unstable. 50 R. Fedorowicz et al. The subsystem consisting of one pipe and two nodes may be solved by a very stable finite difference scheme where the spatial difference of mass density ρ is approximated backwards and the same difference of q – forwards. Calculation of such element can be treated as a block with input and output. The attempt to construct a system of interacting boxes did not give good results, even for the network with a tree-structure. The Gauss-Seidl method applied to the system (1.14) did not give good results. 4.3. The method of solution for networks with tree-structure The method may be applied when the local pressure drop is not taken into account. The finite difference scheme, where the spatial difference ofmass density ρ is approximated backwards and the same difference of q is approxi- mated forwards has been applied. Another link of the unknowns ρ and qwith the unknowns x was applied having in mind that the pressure is continuous at the nodes and the relation between flow ratio in the node enabled reduction of the number of unknowns. This link is especially effective for short pipes. The method has been applied in the previous version of the packet (see [5], Part I). 4.4. The packet development 4.4.1. Forecast module The works are conducted aimed at designing a program or a module of program for loads forecasting. The data concerning loads and the weather are gathered. The base of aggregated and linked data for towns and reduction stations of the first step is prepared. The programwhich enables visualization of the relations between the weather factors, calendar factors and the daily loads Qd and the profiles (see Section 2.4, Part I) is ready ([7], Part I). Theprogramcalculating the load-forecast on thebasis of approximationby hyperplanes forecasts Qd with the error ≈ 6%. The program also calculates the coefficients specific for each day of a week. The preliminary results are more precise for the cold days. The final program will use the neural nets methods. 4.4.2. Improvements • The option, enabling the choice of the first node of the subnetwork by the program operator is rather necessary. For example, whenwe analyze Flow modelling in gas transmission networks... 51 the survival time of the subnetwork, it is very useful to choose the closed valve as the first node. • In the pipes near to the outlet of the compressor station, the energy equ- ation should be taken into considerations. In other pipes, the precision of the isothermal model is quite good. • The formulae for the stationwith reciprocating compressionunits should be implemented. • Thepacket performs calculations when the composition of the gas in the network is constant. The problemof flowwhen the composition of gas in the network changes, should be solved and implemented. The knowledge about the vanishing influence of initial state will be useful in this case. • The calculationswith adaptive time and space step should be implemen- ted. It is important if wewant to automatize the calculations with initial conditions taken from the results of previous calculations and when the fast changes of the flow occurs. • Themodelling of big reservoirs of gas is rather precise. The big reservo- irs are provided with compressor stations. The amount of gas in these reservoirs changes with the seasons of the year, not with days. The big reservoirs are modelled as compressor stations inside the network. The programmodelling phenomena in such gasholders is rather autonomous to the programmodelling the flow in the network. • Modelling thenetworkswith small gasholders is very interesting.A small gasholder without any elements of active control has its state characte- rized by pressure, temperature, mass of the gas and its own parameters and elements. The parameters of flow in different nodes are mutually dependent, so the boundary conditions incorporate different nodes. • The local pressure loss depends on the flow direction in the node. We should take this fact into account when the boundary conditions in two- way nodes are formulated. 4.5. Precise modelling of the station with centrifugal compressors There exists packets modelling precisely, even very complicated stations with centrifugal or reciprocating compressor units (e.g. [2], Part I). The pac- ketmodelling flow in transmission networks should take into account the com- pressor efficiency which depends on thework point P (see Fig.5, Part I). The 52 R. Fedorowicz et al. compressor efficiency influences the temperature at the outlet and the amount of gas used by the turbine driving the compressor unit. When we calculate, successive approximations are necessary. The work-point P depends on the flow and the flow depends on the work-point of the compressor. The calcu- lations performed for a very primitive pipeline with a centrifugal compressor disclosed that the iteration method is divergent, but the method based on the falsi division is effective. The problem is fascinating and we are going to implement the method into the packet. Acknowledgement This work is supported by the State Committee for Scientific Research under the grant No. 0T00A06019. A. Gauss elimination with back substitution and row pivoting A.1. Data int J, double precision R[1 :J], x[1 : J] Class CRROListmembers are: int col, double precision a, pointer to next object of CRROList element. One-dimensional table of pointersCRROList * P[1 :J]. P[i] points at the first element of the i row of data. In the Row-wise Representation Ordered List the elements are ordered up with the col member in each row. A.2. Method for int i=0, step up 1, to i equal to (J−1) • Pivot finding and row pivoting • Elimination and R table modification end of for-loop. A.3. Pivot finding and row pivoting For given i. double precision big= |P[i].a|, int investigate[1 : J−1], int howmany=0, pivot= i; for int r= i, step up 1, to r equal to J Flow modelling in gas transmission networks... 53 if(P[r].col equal to i) 〈 howmany=howmany+1, investigate[howmany] = r if(|P[r].a|>big)〈 big= |P[r].a|, pivot= r〉〉 end of for-loop if(pivot not equal to i) 〈Change P[i] with P[pivot]. Change R[i] with R[pivot].〉. A.4. Elimination and R table modification We know the diagonal identifier i, howmany, and the table investigate. for int h=1, step up 1, to h equal to howmany r = investigate[h]. Execute Arj = Arj −Ari(Aij/Aii) operation for all j > i existing in i-row and r-row. Execute Rr =Rr−Ari(Ri/Aii). end of for-loop. A.5. Conclusions Themost troublesome is the Gauss subtraction Arj =Arj −Ari(Aij/Aii) procedure. Itmust not be against the rule that the element of the resulting list with greater colmembermust be after the element of the list with smaller col member. At each step i the determinant not equal to 0 condition is verified. The procedure does not contain the unnecessary operations of changing the lower triangle of the matrix into zero. The back substitution is presented e.g. in Flannery et al. (1998) for full matrices and it is easy to apply for the sparse matrix in the presented conditions. References 1. Duff I.S., Erisman A.M., Reid J.K., 1986,Direct Methods for Sparse Ma- trices, Clarendon Press, Oxford 2. FedorowiczR., 2001,Gas transportationmodelling inhighpressurenetworks (in Polish), Dissertation,Military University of Technology,Warsaw 3. Fedorowicz R., Kołodziński E., Solarz L., 2002, Flow modelling in gas transmissionnetworks,Part I –Mathematicalmodel,J.Theoretical andApplied Mechanics, 40, 4, 873-894 4. FlanneryB.P., PressW.H., Teukolsky S.A., VetterlingW.T., 1998, Numerical Recipes, The Art of Scientific Computing, Second Edition, Internet, http:/www.nrcom 54 R. Fedorowicz et al. 5. Komarec R., 1999, Transient flow modelling in gas transporting networks with a compressor station (in Polish),Master of Science thesis,WAT (Military University of Technology),Warsaw, June 6. Pissanetzky S., 1984, Sparse Matrix Technology, Academic Press, London Modelowanie przepływu gazu w sieci przesyłowej. Część II – Metoda rozwiązania i algorytmy Streszczenie W części I został sformułowanymatematyczny problem adekwatny do problemu fizycznego, realnie występującego w sieciach przesyłowych gazu. W tej części przed- stawiono rozwiązanie tego problemu. Składa się ono z częścimatematycznej, gdzie dla postulowanego procesu izotermicznego przedstawiono szczegółowe wzory bezwymia- rowe, zaproponowanodyskretyzację przestrzenną i czasową oraz otrzymane równania algebraiczne nieliniowe. Zaproponowano metodę iteracyjną, gdzie na każdym kroku iteracji rozwiązujemy układ równań algebraicznych liniowych wykorzystując fakt, że macierz układu jest rzadka. Pokazano, jak wiedzę o sieci, płynącym gazie, sposobie sterowania i poborach gazu zawrzeć w programie. Szczegółowo pokazano, jak dla dowolnej sieci komputer buduje układ równań algebraicznych, a następnie go rozwiązuje. Mamy nadzieję, że dla badających ten problem publikacja będzie stanowić bazę, do której będą nastę- powały krytyczne odniesienia wskazujące na modyfikację opisu, metod rozwiązania i algorytmów. Rezultaty wskazujące na aplikacyjny charakter pakietu realizującego modelowanie były już publikowane.W kolejnej części przedstawimy uzyskane za po- mocą pakietu rezultaty dotyczące również aspektów fizycznych i obliczeniowych. Manuscript received March 28, 2002; accepted for print July 3, 2002