Keywords: natural gas, gas distribution, pipeline network, flow friction, hydraulic resistance natural gas, gas distribution, pipeline network, flow friction, hydraulic resistance Nomenclature p

Accent is on determination of appropriate friction factor, and on selection of representative equation for natural gas flow under presented conditions in the network.

Calculation of presented looped gas-pipeline network is done according to principles of Hardy Cross method. The final flows were calculated, for known pipes diameters and nodes consumptions while the flow velocities through pipes have to stand below certain values. In optimization problem flows are treated as constant, while the diameters are variables.

When a gas is forced to flow through pipes it expands to a lower pressure and changes its density. Flow-rate, i.e. pressure drop equations for condition in gas distribution networks assumes a constant density of a fluid within the pipes. This assumption applies only to incompressible, i.e. for liquids flows such as in water-distribution systems for municipalities (or any other liquid, like crude oil, etc.). For the small pressure drops in typical gas distribution networks, gas density can be treated as constant, which means that gas can be treated as incompressible fluid [START_REF] Pretorius | A flow network formulation for compressible and incompressible flow[END_REF], but not as liquid flow. Liquid flow and incompressible flow are not synonyms. Under these circumstances flow equation for water or crude oil cannot be literally copied and applied for natural gas flow. This means that original Darcy-Weisbach equation cannot be used without some modifications.

Each pipe is connected to two nodes at its ends. In a pipe network system, pipes are the channels used to convey fluid from one location to another. The physical characteristics of a pipe include the length, inside diameter and roughness. The Darcy's coefficient of hydraulics resistance is associated with the pipe material and age, but also with fluid flow rate and pipe diameter, i.e.

with relative roughness and Reynolds number. When fluid is conveyed through the pipe, hydraulic energy is lost due to the friction between the moving fluid and the stationary pipe surface. This friction loss is a major energy loss in pipe flow and is a function of relative roughness and Reynolds number as mentioned before. In the case when the relative roughness is negligible typical flow regime is hydraulically 'smooth' where Darcy's coefficient of hydraulics resistances depend only on Reynolds number. This regime is typical for gas networks which are built using polyethylene pipes. For that regime Darcy's coefficient of hydraulics resistances can be calculated after so called Blasius type equations or after so called Prandtl type of equations. Prandtl type of equations which is also known as NPK type (Nikuradse-Prandtl-von Karman) is implicit in friction coefficient [START_REF] Moody | Friction factors for pipe flow[END_REF][START_REF] Coelho | Considerations about equations for steady state flow in natural gas pipelines[END_REF]. For liquid flow where the liquid density has larger values compared to gas, Reynolds number increase and if it is accompanied with increased value of relative roughness which is typical for steel pipes, full turbulent regime is most possible. For the fully turbulent regime most convenient is von Karman type of equation from western practice, or Shifrinson equation from Russian practice. For the fully turbulent flow regime Darcy's coefficient of hydraulics resistances depend only on relative roughness.

In this paper is shown new facts in comparison to previous calculation of gas distribution network in Serbian town Kragujevac which was done in 1994. After the implementation, measurements in situ have performed and real measured values deviate from calculated. Previous results are available since they are published and hence comparisons are possible [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF].

This paper addresses to the problem of hydraulic resistance in pipes used for construction of networks for distribution of natural gas in the cities and with subject to all the practical requirements for the engineers charged with design and/or analysis of such system (Mathews and Köhler, 1995). This paper is especially addressed to those engineers willing to understand and interpret the results of calculation properly and to make good engineering decision based on this subject. To predict whether flow will be laminar, hydraulically 'smooth', partially turbulent or fully turbulent, it is necessary to explore the characteristics of flow (Figure 1). Hydraulically 'smooth' regime is also sort of turbulent regime. In these considerations has to be very careful because some of the authors use Darcy's friction factor while the others use Fanning's factor (Brkić 2009a). The Darcy's friction coefficient is four times larger than Fanning's while the physical meaning is equal. Graphically, friction factor for known Reynolds number and relative roughness can be determined using well known Moody diagram [START_REF] Moody | Friction factors for pipe flow[END_REF]. The Darcy friction factor and the Moody friction factor are synonyms. Note also that relative roughness is somewhere defined using pipe diameter and somewhere using pipe radius which can be source of errors in the caseof inappropriate use (Chen,1979;[START_REF] Chen | Comments on 'An explicit equation for friction factor in pipes[END_REF][START_REF] Schorle | Comments on 'An explicit equation for friction factor in pipes[END_REF], Brkić 2009a).

On the Darcy's coefficient of flow friction

Figure 1. Physical description of flow regimes

Note that the Darcy friction factor is defined in theory as λ=(8·τ)/(ρ·v 2 ).

As mentioned in introduction, for polyethylene pipes absolute roughness k is very small compared to the pipe diameter D in i.e. relative roughness is negligible (k/D in → 0), and Darcy's friction coefficient depends only on Reynolds number. For the low value of Reynolds number, but above 2320, flow regime is so called hydraulically 'smooth' (there is no effect of roughness).

For the Reynolds number bellow 2320 regime is laminar. Upper limit for hydraulically 'smooth' regime ζ=16. Typical partially turbulent regime is occurred for 16<ζ<200, and for ζ>200 fully turbulent regime is most possible. Parameter ζ can be found after (1): (1) So, as Reynolds number increases, the flow becomes transitionally rough, i.e. flow regime is partially turbulent in which the friction factor rises above the 'smooth' value and is a function of both, relative roughness and Reynolds number. As Reynolds number increases more and more, the flow reaches a fully turbulent or so called 'rough' regime in which the Darcy's friction coefficient is independent on Reynolds number and depends only on relative roughness. In a hydraulically 'smooth' pipe flow, the viscous sub-layer completely submerges the effect of roughness on the flow. For turbulent flow in 'smooth' pipes friction losses are completely determined by Reynolds number. In rough pipes, however, the value of friction coefficient depends for large values of Reynolds number also on the roughness of the inside pipe surface.

The important point is not so much the absolute roughness size because for the same absolute roughness, the flow resistance in large pipe is considerably smaller than the resistance of a small one.

The Darcy' friction factor for hydraulically 'smooth' regime can be determined after Renouard equation (3):

8 . 1 Re 172 . 0 = λ (2)
This equation belongs to so called Blasius type of equations for hydraulically 'smooth' regime.

For fully turbulent regime from Russian literature [START_REF] Sukharev | Investigation of the hydraulic resistance in polyethylene pipelines[END_REF][START_REF] Nekrasov | Hydraulics for Aeronautical Engineers[END_REF]) is available Shifrinson equation (3):

F o r P e e r R e v i e w O n l y 7 4 1 in D k 11 . 0         ⋅ = λ
(3) Shifrinson equation ( 3) was used for calculation of gas network in Kragujevac in 1994. In our case for Kragujevac gas network, gas dynamic viscosity is η=1.0758•10 -5 Pas which typical for natural gas, density of natural gas is 0.84 kg/m -3 (that implies that relative density is 0.64). Value of absolute roughness is 0.007·10 -3 m for the polyethylene pipes in Kragujevac as reported in the paper of [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF] which is even smaller compared to the value reported in the paper of [START_REF] Sukharev | Investigation of the hydraulic resistance in polyethylene pipelines[END_REF]. [START_REF] Sukharev | Investigation of the hydraulic resistance in polyethylene pipelines[END_REF] found that absolute roughness of inner surface of polyethylene pipes is 0,002·10 -2 m. To find which flow regime are occurred in network it is necessary to find parameter ζ, and hence every tenth pipes from network shown in figure 3 are examined and results are listed in table 1: Table 1

. Determination of flow friction regime

Results of random samples of pipes shown in table 1 concluded assumption that flow in presented gas network is in hydraulically 'smooth' regime because ζ<16, and hence Renouard equation ( 2) is more suitable for this calculation. Shifrinson equation (3) used previously in 1994 which is suitable for 'rough' pipes cannot be used. In some pipes such as 21, regime is not even 'smooth', it is rather laminar. Velocity of gaseous fluids depends on the pressure in pipe since they are compressible (4): Assumption of gas compressibility means that it is compressed and forced to convey through pipes, but inside the pipeline system pressure drop of already compressed gas is small and hence further changes in gas density can be neglected. This is main difference between liquid and incompressible flow. According to this, water flow in pipelines is liquid incompressible flow, while the gas flow is gaseous incompressible flow. Calculated Darcy's friction factor after Shifrinson relation (3) used in 1994 is more than 3 times smaller in comparisons to the results obtained using recommended Renouard equation (2) (Figure 2). [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF] are not in correlation with diagram of [START_REF] Moody | Friction factors for pipe flow[END_REF].

π ⋅ ⋅ = π ⋅ ⋅ ⋅ ⋅ = 2 in 2 in st st D Q 4 D p Q p 4 v (4)
Someone can conclude that bellow laminar and 'smooth' regime, or better to say bellow lines in diagram of Moody which represent these regimes exist some other regimes (Figure 2). Of course, such regimes such as 'sub-laminar' or some kind of turbulent regime below hydraulically 'smooth' regime cannot exist (Figure 1).

General equation for fluid flow

Losses of energy, or head (pressure) losses depend on the shape, size and roughness of a channel, the velocity and viscosity of a fluid, and they do not depend on the absolute pressure of the fluid.

For gaseous fluids some law of thermodynamics also have to be included.

General equation for liquid flow

Experiments show that in many cases pressure drop are approximately proportional to the square of the velocity (5). Equation ( 5) is called the Darcy-Weisbach equation, named after Henry Darcy, a French engineer of the nineteenth century, and Julius Weisbach, a German mining engineer and the scientist of the same era.

ρ ⋅ ⋅ ⋅ λ = - 2 v D L p p 2 in 2 1 (5)
In previous equation velocity and gas density must be correlated, since the gas is incompressible fluid, and hence for gas is more suitable equation in next form (6) because Q·ρ=Q st ·ρ st .

st 2 2 st 5 in 2 2 5 in 2 1 Q 8 D L Q 8 D L p p ρ ⋅ π ⋅ ⋅ ⋅ λ = ρ ⋅ π ⋅ ⋅ ⋅ λ = - (6)
For example pressure drop using Darcy-Weisbach equation for liquid flow (6) in random set of pipes chosen from gas network from Kragujevac is shown in table 2. Table 2. Pressure drops using Darcy-Weisbach equation for liquid flow In (6), if flow rate Q is given for pressure in gas pipeline, i.e. for 4·10 5 abs, and not for normal conditions, density also have to be adjusted for this existing value of pressure in pipeline (volume of gas is four time smaller in 4·10 5 then in 1·10 5 ).

General equation for steady-state flow of gas

Density of gas can be noted as (7). Considering the momentum equation applied to a portion of pipe length, inside which flows a compressible fluid with an average velocity, for example natural gas, and assuming steady state conditions, general equation for gas flow can be written as (8):

0 2 v D L 2 p p T R z M 2 v D dL dp 2 v D dL dp 2 2 in 2 2 2 1 avr avr 2 1 2 2 in 2 1 2 1 2 1 2 in = ρ ∆ λ + - ⋅ ⋅ = ρ λ + ρ = ρ λ + ∫ ∫ ∫ ∫ (8)
In ( 9) flow can be used instead of velocity (8) and combined with (7) gives:

( ) 2 st 2 2 st 2 2 st 2 4 in 2 st 2 st 2 2 st 2 2 2 2 T R z M p D Q 16 A Q A Q v ⋅ ⋅ ⋅ π ⋅ ⋅ = ρ = ρ = ρ (9)
Considering that gas density (see eq 7) at standard pressure conditions is equal as in average pressure in pipeline (ρ st =ρ avr ), and finally assuming that for perfect gas M=M air ·ρ r , general equation for steady-state flow of gas can be written as (10):

st st st air 2 st 2 5 in 2 st 2 2 2 1 T R z M p D Q L 16 p p ⋅ ⋅ ρ ⋅ ⋅ π ⋅ ⋅ ∆ ⋅ λ = - (10)
This equation is rearranged by [START_REF] Renouard | Nouvelle méthode pour le calcul des réseaux maillés de conduites de gaz[END_REF] in very today well known equation (11):

82 . 4 in r 82 . 1 st 2 2 2 1 D L Q 4810 p p C ρ ⋅ ⋅ ⋅ = - = (11)
Previous equation is correlated with Renouard equation used for calculation of Darcy' friction factor for hydraulically 'smooth' regime (2), but with only difference in that the Darcy's factor do not need to be calculated as for equation (10) because it is already incorporated by setting of appropriate coefficient and exponents in equation ( 11). [START_REF] Renouard | Nouvelle méthode pour le calcul des réseaux maillés de conduites de gaz[END_REF] assumed that dynamic viscosity of natural gas is η=1.0757·10 -5 Pas, factor of compressibility is z=1 and that pressure and temperature in pipeline is actually standard temperature (T st =T avg =288.15 K, P st =1.01325·10 5 Pa). This means that by fixing the value of gas kinematic viscosity, the density is also kept fixed, which is physically inaccurate when considering compressible gas flow at medium or high pressure, because the kinematic viscosity of gases is highly dependent upon pressure. According to this, gas flow in city distribution network can be treated as incompressible, like in Kragujevac where pressure drops is minor, and hence Renouard equation ( 11) can be used. Some other transformations of general equation for steady-state flow of gas are available in the paper of [START_REF] Coelho | Considerations about equations for steady state flow in natural gas pipelines[END_REF].

Comparison of the actual and previous results

In figure 3 is shown rings-like part of the gas distribution network from one municipality of the Serbian town of Kragujevac [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF]. There are 29 independent nodes in the ring-like network, 43 branches belonging to rings with 25 branches mutual to the two rings. The total network gas supply input in node 1 is 2339.4 m 3 /h.

Figure 3. Gas-pipeline network in Kragujevac

Hardy Cross method procedure can give good results in design of a looped gas-pipeline network of composite structure [START_REF] Cross | Analysis of flow in networks of conduits or conductors[END_REF][START_REF] Corfield | Distribution design for increased demand[END_REF]Brkić 2009b). Results in table 3 are calculated for locked up diameters from previous calculation available in [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF], while flows are treated as variables. In last two columns in table 3, optimization problem take place where flows are locked up while diameters are variables. Gas consumption per node is presented in figure 3 in brackets and computation results using Renouard's equation adjusted for gas pipelines (11) are shown in Table 3.

. Computational results for gas network in Kragujevac

The flow direction in branches15 and 16 is opposite to the flow direction shown in figure 3.

The velocity limits are 6 m/s for the pipes of small diameter (up to 90 mm) and 12 m/s for the pipes of large diameter (up to 225 mm) according to original project [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF].

But calculated velocities in each of the pipes do not reach even 3 m/s. For example, velocity in some of the pipes is higher than proposed only with assumption that gas is actually a liquid. In previous calculation fact that gas is actually compressed and hence that volume of gas is decreased, is neglected. Hence, mass of gas is constant, but volume is decreased while gas density is according to this, increased. According to this gas network is not nearly optimized and for gas network with input pressure 4·10 5 Pa abs i.e. 3·10 5 Pa gauge all values of velocities are in the ranges of proposed limits. It is much different with liquid. For example gas velocity in pipe 1 from table 1 is only 1.94 m/s, but in the case of water or crude oil equivalent volume of fluid cannot be initially compressed and hence this observed liquid is forced to convey with increased velocity. For pipe 1 from table 1 this velocity for liquid flow is 7.75 m/s.

Optimized design of a gas-distribution pipeline network

In the problem of optimization of pipe diameters, flow rates calculated in previously and shown in table 3 are not any more treated as variable. Results of optimization problem are shown in last two columns in table 3. These flow rates in the next calculation will be locked up, while the pipes diameters will be treated as variable ( 12

L Q 4810 82 . 4 D D L Q 4810 D C D p p ρ ⋅ ⋅ ⋅ ⋅ - = ∂         ρ ⋅ ⋅ ⋅ ∂ = ∂ ∂ = ∂ - ∂ (12)
For example, for the loop XIV from the network shown in figure 3, previous, so called 'loop' equation can be written as ( 13):

( ) { } ( ) ( ) ( ) ( ) { } ( ) ( ) ( ) ( ) ( ) ( )         ⋅ + ⋅ - ⋅ - ⋅ ρ ⋅ ⋅ - = = ∂ ∆ + - - ∂ = ∂ ∂ in 5.
D L Q D L Q D L Q 4810 82 . 4 D D p D C D C D D C (13)
Optimized design of a gas-distribution pipeline network for Kragujevac gas network is shown in last two columns in table 3. This optimization is for average velocity of gas in network of 9 m/s. Some pipes have very small diameters after optimization and hence after adoption of first larger standard diameter, calculation of flow distribution for these new standard diameters has to be repeated.

According to principles of the Hardy Cross method, 'loop' equations, i.e. condition after second Kirchhoff's law has to be fulfilled at the end of calculation. These equations represent energy continuity, while the first Kirchhoff's law represents mass continuity for nodes. Mass continuity has to be fulfilled in all iterations for all nodes without exceptions. Diameter corrections are calculated after (14):

{ } ( ) ( ) { } i i i i D D C D C D ∂ ∂ = ∆ (14)
These corrections calculated for all particular loops from the Kragujevac gas network have to be added to previous value of diameter according to algebraic scheme available in literature [START_REF] Corfield | Distribution design for increased demand[END_REF]Brkić 2009b). Of course, many more efficient methods than original Cross (1936) method exist, but this is not main subject of this article. For example, calculated flows presented in fourth column in table 3 are obtained using the node-loop method [START_REF] Boulos | Comprehensive Water Distribution Systems Analysis Handbook for Engineers and Planners[END_REF] and result of optimization problem from the last two columns in table 4 are obtained using improved Hardy Cross method [START_REF] Boulos | Comprehensive Water Distribution Systems Analysis Handbook for Engineers and Planners[END_REF]Brkić, 2009b). In both mentioned method matrix formulations are used which can be easily solved using MS Excel.

Conclusions

Hardy Cross method procedure or similar improved sort of procedures can give good results in design of a looped gas-pipeline network of composite structure. According to the price and velocity limits, the optimal design can be predicted. But all parameters, such as friction factor, relation for calculation of pressure drop in pipes i.e. equation for calculation of gas flow must be chosen in a very careful way. Today, distributive gas network is usually calculated using Renouard equation for determination of gas flow and pressure drop value. Inner surface of polyethylene pipes which are almost always used in gas distribution networks are practically smooth and hence flow regime in the typical network is hydraulically 'smooth'. Using of inappropriate friction factor can lead to a paradoxical result, i.e. calculated Darcy's friction factor can belong to nonexistent regimes such as 'sub-laminar' or turbulent regime which is bellow hydraulically 'smooth'. This means that calculated Darcy's friction factor are highly underestimated. Further using Darcy-Weisbach equation for liquid flow instead it modified version for gaseous flow, pressure drop are overestimated. This leads to pseudo-accurate final results. After all, visible error is occurred in numerical values of some parameters, like velocity in this case. Velocities in network according to previous results are significantly larger than Hazen-Williams relation is frequently used for waterworks or sewerage systems [START_REF] Boulos | Comprehensive Water Distribution Systems Analysis Handbook for Engineers and Planners[END_REF] but not only inaccurate, or better to say accurate in limit range, the Hazen-Williams equation is conceptually incorrect [START_REF] Liou | Limitation and proper use of the Hazen-Williams equation[END_REF]. Even assumed that gas flow is incompressible in municipality pipelines, this relation cannot be used for such systems as gas-distribution networks. Table 1. Determination of flow friction regime

Figure 2 .

 2 Figure 2. Identification of the problem in a characteristic diagram

  is that diameters of pipes are too large for required amount of gas and hence network is not optimal.In the project for gas pipeline in Serbian town Kragujevac from 1994, the Darcy-Weisbach equation for liquid flow were used instead of modified Darcy-Weisbach equation for gaseous flow, accompanied with inappropriate usage of Shifrinson equation for full turbulent regime which is typical for liquid flow in a steel pipes, instead of some sort of Blasius or Prandtl form of equations which are typical for flow of gas through polyethylene pipes which are used for Kragujevac gas network.

  Figures

Table 2 .

 2 Pressure drops using Darcy-Weisbach equation for liquid flow

Table 3 .

 3 Computational results for gas network in Kragujevac

Table 1 .

 1 Determination of flow friction regime Renoard equation can be used if ζ<16, and Shifrinson if ζ>200, but in both cases only if Re>2320

	Pipe	D in (mm)	Q (m 3 /s)	Velocity (m/s) liquid a gas	Re	Friction coefficient λ b Renouard c Schifrinson	d Criterion -ζ Renouard Schifrinson
	1	220.4	0.295497222	7.75	1.94	20391.6 0.028828830	0.008258	0.11	0.06
	11	158.6	0.126230556	6.39	1.60	12105.2 0.031666022	0.008966	0.10	0.05
	21	15.4	0.000383333	2.06	0.51	378.6	0.05908245	0.016062	0.04	0.02
	31	96.8	0.040133333	5.45	1.36	6305.8	0.035610262	0.010144	0.09	0.05
	41	96.8	0.059622222	8.10	2.03	9367.9	0.033161384	0.010144	0.12	0.07
	51	109.8	0.086647222	9.15	2.29	12002.3 0.031714738	0.009829	0.14	0.08
	61	79.2	0.021661111	4.40	1.10	4159.7	0.038379256	0.010666	0.07	0.04
	a for 4 bar abs; see (4)
	b see (2)
	c see (3)
	d see (1);

Table 2 .

 2 Pressure drops using Darcy-Weisbach equation for liquid flow

	Pipe	L (m)	a Pressurre drop Pa Renouard Schifrinson
	1	84	276.84	79.30
	11	119	407.40	115.35
	21	212	1446.82	393.32
	31	115	528.42	150.52
	41	278	2625.35	803.07
	51	78	792.36	245.57
	61	383	1506.96	418.78
	a using (5) or (6) and values from Table 1; gas density
	is ~0.84

Table 3 .

 3 Computational results for gas network in Kragujevac

	Optimized design

a first larger or smaller standard diameter have to be chosen

URL: http://mc.manuscriptcentral.com/lpet Email: jamessp8@aol.com Petroleum Science and Technology


Towards multiscale archival digital data 1

Laurent Romary (Inria & DARIAH), Charles Riondet (Inria)

In this paper, we would like to present some ideas on the use of the archival standards in various contexts that exemplify the complexity of such standards and provide users with innovative ways to handle EAD content. Our main idea is that researchers, Cultural heritage institutions, archival portals and standards maintenance bodies could greatly benefit from a multiscale modelling of archival data, but also from multiscale representations and documentations. A first step is on the way to being cleared in the domain of the management of heterogeneous archival sources in one single environment, namely a federated portal, like in EHRI. We built a methodology based on a specification and customisation method inspired from the long lasting experience of the Text Encoding Initiative (TEI) community. In the TEI framework, one has the possibility of defining project-specific subsets or extensions of the TEI guidelines while maintaining both the technical (XML schemas) and editorial (documentation) specification within a single framework. Using the same framework for EAD data allows us to express precise content-oriented rules combined with some interesting possibilities of integrating the human readable documentation in the validation process.

EAD, What's wrong with it

The development of EAD was initiated in 1993 at the Library of Berkeley, with the idea of building a non proprietary format for finding aids, reflecting the hierarchical structuration of archival fonds. If preliminary attempts were expressed in SGML, the first version of EAD used XML, and was released in 1998. A second version was released soon after in 2002, EAD2002, which is still the most widely used version. It is maintained by the Library of Congress and the Society of American Archivists. In 2010, a global revision process was initiated, in order to make EAD more connected to Linked Data technologies, and to reach a better integration with the others XML archival formats: EAC-CPF and EAG and in 2015, EAD3 was officially released. However, in the world of cultural heritage institutions and research, archival description is often considered as a pending issue, a hindrance to data exchange and accuracy. Since its creation, EAD faces criticism, as many observers are pointing to its permissiveness as a problem. Yet in 2001, Shaw asks for a "more prescriptive descriptive standard" [START_REF] Shaw | Rethinking Balancing Flexibility and Interoperability[END_REF].

Still today, and even if EAD32 is globally seen as a step in the right direction, EAD is generally seen as a poorly structured and interoperable standard, not very suitable for data exchange, and is paradoxically considered by some information specialists, not a "standard for archival description" [START_REF] Bunn | Developing Descriptive Standards: A Renewed Call to Action[END_REF]. We will not go any further in this controversy, but point the fact that the archival community, though aware of these weaknesses, still broadly works with EAD and is still willing to improve the quality of digital archival descriptions.

There is room to improve EAD in two main aspects: 1) handle its flexibility and 2) preserve all the complexity of the content when exchanging archival description. Of course, the new Records in Context content mode [EGAD 2106] proposes a nice way to handle these issues, with an ontology meant to bring together all the pieces of archival information (authorities, institutions, functions and records), natively compliant with semantic web technologies. But, before this solution is adopted and implemented, EAD still is and will be the archival community standard. The framework we propose will allow for better exchange and dialog between archival data and together with others resources available online.

The EHRI use case

The EHRI environment is a perfect use case to apply our method, because of the heterogeneity of the corpus, characterized by a great diversity of languages, description levels, and archival practices, and the goal to ingest all these archival descriptions in one single environment. These various sources need therefore to be compared, checked in quality, and processed before being integrated in the repository. To do so, the pivot format is naturally EAD (version 2002), used for automatic ingestion in EHRI database and also for exports. Like for all the archival portals, the two crucial questions are how to deal with so many differents ways of encoding EAD, and how to guarantee that the descriptions are compliant with EHRI requirements. To handle this situation, we propose a method to create customizations for EAD in order to refine archival descriptions both in the structure and in the content, and of course respect entirely the EAD syntax.

This method is developed in the context of the umbrella project Parthenos3 which aims, among other things, at disseminating information and resources about methodological and technical standards in the humanities. One of the main objectives of Parthenos is to create a Standardization survival Kit (SSK) [Romary et al., 2016], whose main features are to:

-Propose generic research scenarios to scholars where the use of standards play a key role -Communicate around community initiatives -Support standardization activities in domains where it is needed.

Within Parthenos, one of the scenarios we will provide in the SSK is precisely a scenario guiding scholars and cultural heritage information specialists in the creation of projectspecific EAD schemas.

The TEI-ODD specification framework

In this project, we are inspired by another very strong community standard : the Text Encoding initiative. This format facilitate the representation of any textual resource in XML. It was built for digital editions of historical texts, but can be used in many other situations. For instance, what we are interested in is a subset of the TEI meant to create XML formats specification (the TEI itself is described with this subset of TEI). This is called "One document does it all" and it allows us to model specific subsets, extensions or profiles of the described format. ODD can be used to refine the behaviour of elements and attributes, for any XML format, contains all the human readable documentation and can be processed to generate various resources: a validation schema (in many formats) and some documentation (in many formats).

ODD is based on the principles of literate programming, which means that this language combines formal (specifications) and informal declarations (descriptive prose and examples) [Knuth 1983]. It combines in the same environment the technical specifications and the user guidelines for the key components of the TEI Abstract Model, primarily elements and attributes, but also modules, classes and macros [START_REF] Lou | RelaxNG with Son of ODD[END_REF]. For example, to write the specification of an element, the tag used is <tei:elementSpec>. It contains elements for documentation, like the <tei:gloss> (a phrase or word used to provide a gloss or definition) or <tei:desc> (a brief description of the object documented by its parent element, typically a documentation element or an entity). The <tei:classes> element is used here to link elements with their attributes, and the <tei:content> contains the relaxNG specification, i.e. what elements can be children of the described element. The official EAD schema and the official EAD tag library were encoded in an ODD document (Agreement of the Library of Congress and the Society of American Archivists), in the context of the Parthenos project4 . This EAD ODD is a starting point for EHRI, used to create an EHRI-specific EAD profile with very precise content oriented rules based on EHRI requirements and on the CHI data models and some qualitative documentation to be served to the user of conversion and validations services provided by the EHRI project.

EHRI specific rules

EHRI has its own ODD, project specific, that inherits everything from the generic EAD ODD, except the elements and attributes that have a different behaviour in EHRI. The philosophy is to keep the EAD schema as it is, i.e. not modify directly the RelaxNG specifications. Instead, we use another validation language: ISO Schematron5 . EHRI already used schematron rules to control the input descriptions. We completed them, respecting the same organisation. Schematron validation serves diagnostics to the content providers, by emphasizing:

-technical errors and proposes a solution, as EHRI conceive it -EHRI descriptions guidelines requirements -EHRI descriptions guidelines proposals, or "nice to have" points Some rules reflects the requirements of EHRI database content model. For instance, it asks that the <date> elements contains a @normal attribute whose content respect the ISO8601 standard on representation of dates and time. This constraint is expressed in the ODD file with embedded schematron in the following way:

<elementSpec ident="date" module="EAD" mode="change"> <constraintSpec ident="dateNormal" scheme="isoschematron" type="EHRI" mode="add"> <desc>All the <gi>date</gi> elements MUST have a <att>normal</att> attribute whose pattern respects the ISO8601 standard and take the following form: YYYY-MM-DD</desc> <constraint> <sch:rule context="date"> <sch:assert role="MUST" test="matches(@normal,'^(

([0-9]|[1-9][0-9]|[1-9][0-9]{2}|[1-9][0-9]{3}))-(0[1-9]|1[012])- (0[1-9]|[12][0-9]|3[01])$')">@normal attribute MUST respect ISO8601 pattern = YYYY-MM- DD</sch:assert> </sch:rule> </constraint> </constraintSpec> </elementSpec>
This second rule is also a requirement, but for different reasons. For the sake of comprehension of the archival description, EHRI requires that a <scopecontent> element should be present somewhere. the choice is let to the provider to write on general paragraph and put it at the highest level (<archdesc>) or add a more precise <scopecontent> for each subcomponents, from <c01> to <c06>. Here, the rule is called at the <archdesc> level, because it is more likely that the CHI provides a global <scopecontent> if it didn't exist before.

<elementSpec ident="archdesc" mode="change"> <!--… --> <constraintSpec ident="scopecontentInArchdescOrC" scheme="isoschematron" type="EHRI"> <desc>A <gi>scopecontent</gi> element SHOULD be present in the description at least in <gi>archdesc</gi>, if not in the <gi>c</gi> elements.</desc> <constraint> <sch:rule context="archdesc" role="SHOULD"> <sch:assert test="scopecontent or dsc/c01/descendant-orself::scopecontent">a "scopecontent" element SHOULD be present at least in "archdesc" if not in the "c" elements</sch:assert> </sch:rule> </constraint> </constraintSpec> </elementSpec>

The last rule showed is the lowest level of constraint. It presents some possibilities to make the description more complete. In particular, these rules focuses on the content related elements of <archdesc>. Therefore, these messages are not considered as real errors, but as pieces of advice that the providers can follow or not.

<elementSpec ident="archdesc" mode="change"> <!--… --> <constraintSpec ident="bibliographyPossible" scheme="isoschematron" type="EHRI"> <desc>The <gi>archdesc</gi> element COULD contain a <gi>bibliography</gi> element.</desc> <constraint> <sch:rule context="archdesc"> <sch:assert role="COULD" test="bibliography">archdesc COULD have a bibliography</sch:assert></sch:rule> </constraint> </constraintSpec>

The rules added to the EAD schema in EHRI specify all the different parts of the archival description : the administrative metadata (the <eadheader>, in particular the history of the modification of the EAD), the description itself (<archdesc>, <c> and <did>), and the content elements (the access points, i.e. the named entities, persons, places, organisations, …, but also the dates). Another type of specific rules is related to all the standardized codes used to identify some pieces of information, like the languages used (ISO639), the archives (ISO15511), etc.

Overview of EHRI schematron rules: The <processinfo> element should contain a <date> element as descendant. SHOULD

A <scopecontent> element should be present in the description at least in <archdesc>, if not in the sub components elements.

SHOULD

The sub components elements should be numbered components between <c01> and <c06>

SHOULD

The <archdesc> element could contain these content related elements:

-<langmaterial> -<custodhist> -<otherfindaid> -<originalsloc> -<altformavail> -<bibliography> -<odd> -<note> -<controlaccess> COULD <langmaterial> could contain a <language> element. COULD
If the element <altformavail> is not empty, you could try to identify if the originals are present in the EHRI portal and make a link between the two descriptions.

COULD

If the element <originalsloc> is not empty, you could try to identify if copies are present in the EHRI portal and make a link between the two descriptions. COULD 3. Entities a. Access points

EHRI Rules Role

In <controlaccess>, EHRI welcomes any access points types : <subject>, <geogname>, <persname>, <orgname>.

COULD

Access points could be chosen in authority lists. The list is declared with a @source attribute. The related identifier of this authority should be declared in an @authfilenumber attribute. Note that EHRI provides URLs for vocabularies and authorities.

COULD

In the access points, person names should be structured like this : Family name, given name SHOULD b. <unitdate> and <date>

EHRI Rules Role

<unitdate> should have a non-empty @normal attribute SHOULD

The @normal attribute of <unitdate> must respect the ISO8601 pattern: YYYY-MM-DD MUST <unitdates> could have a @label attribute or an @encodinganalog attribute, describing the type of date COULD All the <date> elements must have a @normal attribute whose pattern respects the ISO8601 standard and take the following form: YYYY-MM-DD MUST 4. Languages and other coded values

EHRI Rules Role

<language> must have a @langcode attribute, taken from the list given by the ISO639 standard.

MUST <language> must have a @scriptcode attribute, taken from the list given by the SHOULD ISO15924 standard.

If the language of the description is not English, a parallel form of the title in English should be added. For instance, using another <unittitle> element with a @type attribute SHOULD <eadid> should contain a @mainagencycode attribute, which provides (if applicable) the ISO15511 code for the institution that maintains the finding aid.

SHOULD

If the @repositoryencoding is set to iso15511, the format of the value of the @repositorycode attribute is constrained according to the International Standard Identifier for Libraries and Related Organizations (ISIL: ISO 15511): a prefix, a dash, and an identifier.

SHOULD

In the EHRI mapping and validation workflow, the EHRI-EAD schema is used to test the archival descriptions before they are ingested in the portal. The result is of this validation is a list of messages (presented above) linked to precise fragments of the tested description. Therefore, the archive that ingests its descriptions in EHRI portal is informed of the changes it has to make to be sure its data could be integrated in the portal harmlessly. In the future, it is also planned that some uncritical modifications could be automatically made inside the validation framework (base on the Schematron Quickfix extension7 ).

Conclusion

Offering a standard-based method to gain interoperability between heterogeneous data allows users, above all researchers, to access high quality standardized data. On the other hand, a small CHI sharing easily its data via the EHRI portal gains visibility, by showing easily underexposed data, and creates data enrichments opportunities.

This method may be of a wider interest within similar environments (i.e. archives portals). As it is one of the components of the Parthenos Standardization Survival Kit -a solution that offers researchers needing standardized methods and resources complete frameworks to carry out their project, in Arts and Humanities and Heritage science, it can be used freely by any interested project. Parthenos is also willing to give support and maintain the EAD ODD for a substantial period. More, this solution can be seen as a possible bridge between EAD2002 and EAD3, and more broadly could be considered as a tool for the future maintenance of the EAD standard, in order to, like for the TEI, oriente this maintenance towards a (wise) ever ongoing revision methodology.

It could also be an opportunity to bring together EAD and TEI and propose on the fly generation of skeletal TEI documents based on EAD descriptions.

Fig. 1 :

 1 Fig. 1: Example of the ODD specification for an EAD element

  <eadheader> must contain a <profiledesc> element MUST the <eadid> element must contain text. Most of the time, it is automatically generated by the archival tool. MUST <eadheader> must contain information on the language used in the EAD document with a <langusage> element containing a <language> element MUST <eadheader> should contain a <creation> element SHOULD <eadheader> should specify a <publisher> SHOULD The <titlestmt> element should contain a <author> element SHOULD A date of creation for the finding aid is welcome. The relevant element is <date>, child of <creation> COULD <descrules> has a default value added automatically by EHRI. Therefore, the content of <descrules> will be overwritten In <revisiondesc>, each <change> element should contain a <date> element and a <item> element. [This rule has been taken from Library of Harvard Archivesspace Checker 6 ] SHOULD The <date> element for each <change> in <revisiondesc> should not be empty SHOULD 2. Archival description (<archdesc>, <did>, <c>)EHRI Rules Role<archdesc> must have a level attribute. MUSTThe value of the <archdesc> @level attribute should be limited to four values: must have a @type attribute MUST if <dsc>'s @type attribute has "othertype" for value, <dsc> must have a not empty @othertype attributeMUSTThe sub components elements (<c01> to <c06>) must have a @level attribute. MUSTIf the @level attribute of <c01>, <c02>, … has the value 'otherlevel', an attribute @otherlevel MUST be added MUST <did> elements must contain: -a <unitid> -at least on non-empty <unittitle> MUST Each unit of description should have an identifier in the element <unitid>. SHOULD In a given EAD document, all the <unitid> elements must be unique MUST In the <did> element, <physdesc> should come with a non-empty <extent> element SHOULD <archdesc> should contain a non-empty <origination> element. SHOULD <archdesc> should contain a non-empty <processinfo> element. SHOULD

This work is developed in the context of the H2020 projects EHRI and PARTHENOS

https://www.loc.gov/ead/EAD3taglib/index.html

http://parthenos-project.eu

http://github.com/ParthenosWP4/standardsLibrary/blob/master/archivalDescription/EAD/odd/EADSpe c.xml

http://schematron.com/

https://github.com/harvard-library/archivesspacechecker/blob/master/schematron/archivesspace_checker_sch.xml

http://schematron-quickfix.github.io/sqf/publishing-snapshots/April2015Draft/spec/SQFSpec.html


Dejan Brkic
email: dejanbrkic0611@gmail.com

The authors of the discussed paper show a possible strategy for dealing with zero-flows in solving the nonlinear equations for water distribution systems when the Hazen-Williams equation is used. Recently, [START_REF] Elhay | Dealing with zero flows in solving the nonlinear equations for water distribution systems[END_REF] presented a similar method for solution of the zero-flow problem also when the Hazen-Williams model is used, but they also explain and give a solution for the possible problem with zero flow when the Darcy-Weisbach model is used. In this discussion, a few simple remarks how to avoid the zero-flow problem in a network of pipes will be highlighted. Also, possible physical interpretation related to the problem will be explained.

Zero-flow in Hazen-Williams model

Both contributions, by the authors of the discussed paper and by [START_REF] Elhay | Dealing with zero flows in solving the nonlinear equations for water distribution systems[END_REF], to the solution of the zero-flow problem when the Hazen-Williams model is used, cannot be disputed. Mathematical interpretation of the problem from both papers stands, but at the same time everybody has to be aware that the Hazen-Williams equation, used in both papers is obsolete and hence should not be used [START_REF] Liou | Limitation and proper use of the Hazen-Williams equation[END_REF], Brkić 2012a[START_REF] Simpson | Closure to 'Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model[END_REF].

Zero-flow can occur when the Hazen-Williams formula is used since the coefficient is always The fact that the Hazen-Williams model is used for calculation in EPANET is also avoidable since this software equally allows the use of the Darcy-Weisbach model (Simpson andElhay 2011, Brkić 2012a). Because the Darcy-Weisbach model with the Colebrook formula for the friction factor is theoretically more sound (Brkić 2011a[START_REF] Brkić | Determining friction factors in turbulent pipe flow[END_REF], the usage of the Hazen-Williams equation is strongly discouraged. Finally, the Darcy-Weisbach model can be used also for calculation of gas distribution networks, while the Hazen-Williams model cannot in any circumstances [START_REF] Brkić | An improvement of Hardy Cross method applied on looped spatial natural gas distribution networks[END_REF]2011b,c).

Zero-flow in Darcy-Weisbach model

On the other hand, the zero-flow problem can occur when the Darcy-Weisbach formula is used only if laminar flow takes place [START_REF] Elhay | Dealing with zero flows in solving the nonlinear equations for water distribution systems[END_REF], Simpson and Elhay 2011, Brkić 2012a). This is because the resistance is independent of flow when the Darcy-Weisbach formula is in use only in the case of a laminar flow regime. So, knowing that laminar flow can occur only rarely and only in a few pipes of a water distribution network, calculation for these pipes should be perform as for the other pipes in which turbulent flow takes place. Further calculation with this assumption will not introduce significant error in the final result.

Existence of pipes with laminar flow only means that the model of the network is not rationally planned. This subsequently means that diameters of these pipes have to be changed.

Note that the network should be calculated for maximum possible nodal demands, which means that the network is rationally planned only if turbulent flow takes place in all pipes.

Analogy with electrical networks

It is true that laminar flow resistance in the Darcy-Weisbach interpretation is a constant for a single pipe [START_REF] Elhay | Dealing with zero flows in solving the nonlinear equations for water distribution systems[END_REF][START_REF] Simpson | Closure to 'Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model[END_REF], Brkić 2012a). This means that flow resistance, r≠r(λ), in the laminar regime does not depend on the value of the Darcy friction factor, λ (for the laminar regime, the Darcy friction factor can be calculated as λ=64/R, where R is the dimensionless Reynolds number). On the other hand, in the turbulent regime, flow resistance does depend on the Darcy friction factor, i.e. r=r(λ) (where the Darcy friction factor can be calculated using the well known Colebrook formula). To make a point, a clear analogy with electrical resistance exists in the case of resistance in laminar flow. So, knowing that electrical networks can be solved in a non-iterative procedure using only Ohm's and two Kirchhoff's laws, it can be concluded that hydraulic networks can be equally solved using some sort of Ohm's law rearranged for use in hydraulic networks and two Kirchhoff's laws. Laminar flow resistance is independent of flow, but the whole calculation will be spoiled if even a single pipe of the hydraulic network has turbulent flow (a single pipe with turbulent flow renders impossible a non-iterative calculation of the whole network). In such a network, in which in all pipes laminar flow takes place, pipes with zero flow will be treated simply as a break in the circuit (a connection with infinity large resistance) or as a totally choked pipe, which will not cause any problem since no iterative procedure is needed.

Division by "zero" in computer environment

Computers today use the IEEE standard for arithmetic precision and therefore small numbers bellow a standard boundary will also be treated in the computer as zero which also can lead to the singularity of matrices used in calculation of water distribution network (Brkić 2012c, Sonnad and[START_REF] Sonnad | Constraints for using Lambert W function-based explicit Colebrook-White equation[END_REF]. Also, use of software specialized only for matrix calculation (such as MatLab by MathWorks or even MS Excel) can be sometimes recommended as a better solution compared with the use of specially developed software for a water distribution network. In MatLab, it is possible to devise all parts of the calculation, while in a specialized software program for water networks, such EPANET, the designer is more restricted since the calculation procedures are already incorporated in the program code.

Possible physical interpretation of "zero-flow"

Although pipes with no flow in a real looped network of pipe can exist, it is more likely that a quite unrealistic model of a water distributive network is chosen if zero flow occurs (or the model does not accurately represent the system). Considering the network model from Figure 1 of this discussion which has a vertical axis of symmetry (symmetry in pipes diameters and nodal demands). Obviously such a network is excellent for the examination of the zero-flow problem. Symmetric networks can be found in [START_REF] Elhay | Dealing with zero flows in solving the nonlinear equations for water distribution systems[END_REF] and in Álvarez et al. (2011). A symmetric network was referred to in the discussed paper in the work of [START_REF] Elhay | Dealing with zero flows in solving the nonlinear equations for water distribution systems[END_REF]. To further illustrate the point of the shown zero-flow problem, the non-zero demand of node 2 of the network from Figure 1 is equal to with the demand of node 3, node 4 equal to node 5 and node 6 equal to node 7. Also, it can be assumed that all pipes have the same diameter. In that way symmetry of the network and symmetry of node demands leads to the logical conclusion that zero-flow takes place in pipes 2, 6 and 9. This subsequently leads to the conclusion that the consumer connected to pipes 2, 6 and 9 will suffer of water shortage since water users are really located between junctions (Figure 2). In reality, the consumers connected to pipes 2, 6 and 9 will almost certain have enough water since these pipes are supplied from two sides (two-way supplied pipes). Or in other words, the lowest pressure of water is somewhere between the two nodes [START_REF] Brkić | An improvement of Hardy Cross method applied on looped spatial natural gas distribution networks[END_REF]). This situation is not allowed and cannot be calculated using any of the Hardy Cross type methods of s for calculation of looped pipe networks [START_REF] Brkić | Iterative methods for looped network pipeline calculation[END_REF]. For example, the normal situation for pipe 5 is that water flows from node 3 towards node 5. This means that the pressure in node 3 is higher than the pressure in node 5 with a monotonically decreasing pressure through pipe 5.

On the other hand, the pressures in nodes 2 and 3 of the network from Figure 1, are equalized which means that flow through pipe 2 is logically impossible. This assumption can be disputed knowing that the point of the lowest pressure (lower than in nodes 2 and 3) in reality is somewhere between these two nodes. This situation produces simultaneous flow from node 2 towards node 3 and also from node 3 to node 2 (two-way flow or simultaneous flow from two opposite directions). This is possible if the nodes in a model of the network are poorly spatially distributed. A good engineer should know that the real consumers are not concentrated in a node (Figure 2). They are actually distributed between nodes. Consumption concentrated in a node is only a model of the real situation. Also, nodes are not necessarily the only junctions in a network (Figure 3). In the network from Figure 1, nodes should be placed also between nodes 2 and 3, between nodes 4 and 5, and also between nodes 6 and 7 (nodes 9, 10 and 11 in figure 3 of this discussion). The actual situation of the demand pattern will in that way be modeled more realistically (Figure 3). It also has to be noted that an initially poorly conditioned network has as the consequence a poorly conditioned Jacobian which leads directly to a singularity in the related matrix. The general recommendation is that the symmetry in a network should be avoided and if the symmetry exists, nodes at least should be always placed at the axis of symmetry (in that case a node should be placed at every point where pipes and the axis of symmetry cross each other). Symmetry of node demands and pipe diameters also should be avoided.

To conclude, temporary zero-flow rarely can occur in some of the pipes during the calculation of a looped network (virtual change of flow direction during the iterative procedure usually does not cause the zero-flow problem). But, if zero-flow remains as is at the end of the calculation, this usually means that the modeled network is not a good image of the real situation in the field.

  Vol. 139, No. 4, pp. 456-459. DOI: 10.1061/(ASCE)HY.1943-7900.0000694. Dejan Brkić, Ph.D., Research Associate, University of Novi Sad, Faculty of Technology, Bulevar cara Lazara 1, Novi Sad, Serbia, and Strumička 88, 11050 Belgrade, Serbia, e-mail: dejanbrkic0611@gmail.com

  Manuscript DBClick here to download Manuscript: Discussion DB.doc independent of flow. The argument that the Hazen-Williams model can be used since it has been in common use for a very long time[START_REF] Simpson | Closure to 'Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model[END_REF], simply does not stand.

Figure 1 .

 1 Figure 1. Unrealistic symmetric model of water distribution network (chosen only for the

Figure 2 .

 2 Figure 2. Modeled versus possible real situation with two-way flow in a water distributive

Figure 3 .

 3 Figure 3. Good conditioned node pattern in the water distributive network

Figure 1 DB

 1 Figure 1 DB Click here to download high resolution image

Figure 2 DB

 2 Figure 2 DB Click here to download high resolution image

Figure 3 DB*

 3 Figure 3 DB Click here to download high resolution image


Ph.D Dejan Brkić

Discussion of "Jacobian Matrix for Solving Water Distribution System Equations with the Darcy-Weisbach Head-Loss Model" by Angus Simpson and Sylvan Elhay

June 2011, Vol. 137,No. 6,.1061/(ASCE) HY.1943-7900.0000341 Dejan Brkić, Ph.D., Petroleum and Natural Gas Engineer, Strumička 88, 11050 Belgrade, Serbia, e-mail: dejanrgf@tesla.rcub.bg.ac.rs

The general recommendation by [START_REF] Simpson | Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model[END_REF] that the Hazen-Williams formula be replaced by the Darcy-Weisbach is heartily applauded. Today, this fact cannot be disputed.

The main purpose of this discussion is to recommend more appropriate formulas for the calculation of the Darcy friction factor, λ which are more suitable for computation of the derivatives in a Jacobian matrix, and also to clarify a few points in the calculation of the flowdependent Jacobian. The Darcy-Weisbach formula contains the Darcy friction factor, λ which in general is dependent on the Reynolds number and, therefore, in general on flow, i.e. λ(Re)→λ(Q). This fact can be incorporated in the Jacobian matrix [Eq. (6) in the discussed paper] as clearly shown in [START_REF] Simpson | Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model[END_REF]. On the other hand, flow, Q, in the formula for flow friction, λ, is not treated as a variable in EPANET (Computer software, Washington, DC, U.S. Environmental Protection Agency) during the computation of the Jacobian matrix (probably following the analogy with the Hazen-Williams scheme). Use of the Darcy friction factor dependent on flow λ=λ(Q) in the Jacobian matrix is more accurate but also more complicated. A certain balance between accuracy and complexity in engineering systems has to be carefully evaluated. The procedure shown in the discussed paper is more accurate mathematically and also can be more easily implemented in a computer code. This increased computing burden of the improved procedure has to be Text DB Click here to download Manuscript: Discussion DB (Simpson and Elhay).doc carefully investigated where the ratio between complexity and accuracy will decide whether this method will be used in every day engineering practice [START_REF] Giustolisi | Some explicit formulations of Colebrook-White friction factor considering accuracy vs. computational speed[END_REF].

The pressure drop model Δp (Pa) for the Darcy-Weisbach formula Δp=8

•L•ρ•Q 2 •λ/(π 2 •D 5 ) is
used in this discussion instead of the head-loss model ΔH (m) which was used in the discussed paper where the connection between them is simple, i.e. Δp=ρ•g•ΔH. Also, the Darcy friction factor noted as λ is used as in [START_REF] Colebrook | Turbulent flow in pipes with particular reference to the transition region between the smooth and rough pipe laws[END_REF] where the Darcy friction is noted as f in the discussed paper (one should not confuse the Darcy friction factor used in this discussion and the discussed paper with the Fanning friction factor which is numerically equal to one-fourth the Darcy friction). The Darcy friction factor λ for the turbulent zone (turbulent flow is the most common in pipe networks) can be calculated using the implicit Colebrook Eq. ( 16), using some of the explicit approximations to the Colebrook equation (Brkić 2011a) or using some other suitable formulas [EPANET uses the approximation by [START_REF] Swamee | Explicit equations for pipe flow problems[END_REF] which is shown in the next formula on the far right side of the equation].

                       D D 7 . 3 Re 74 . 5 log 2 7 . 3 Re 51 . 2 log 2 1 9 . 0     (16)
The first derivative of the Darcy-Weisbach equation (appropriate terms in the Jacobian matrix [Eq. ( 6) in the discussed paper]) has to be calculated using the fact that the Darcy friction factor, λ contains flow, Q, which now has to be treated as a variable, i.e., λ=λ(Q) and not as constant anymore (e.g. flow is treated as constant in [START_REF] Brkić | Spreadsheet-based pipe networks analysis for teaching and learning purpose[END_REF] and in EPANET). Explicit

    Q D D Q D Q L Q D Q L Q Q p                                                                           2 9 . 0 10 5 2 5 2 7 . 3 4 74 . 5 log 16 16         (17)
This procedure is not simple and hence the practical achievement in accuracy can be disputed knowing the empirical nature of the Colebrook and similarly-related formulas (with an error of a few percent or more compared with experimental data).

Further, the Hagen-Poiseuille equation used for the laminar regime in the discussed paper [Table 2 and Eq. (2) in the original paper] is nothing more than the Darcy-Weisbach formula where the suitable Darcy friction factor λ=64/Re is included (Eq. 18).

Q r D Q L D Q L p                4 5 2 2 128 Re 64 8      (18)
The fact that hydraulic resistance, r, from Eq. ( 18) is independent of flow [as shown in Eq. (2) of the discussed paper] implies that a water network in which laminar flow takes place (in all pipes) can be solved using a non-iterative procedure. On the other hand, the Darcy friction factor, λ in the laminar zone depends only on the Reynolds number, Re, which further depends on flow Q. The fact that laminar resistance, r is independent of flow, Q, is not crucial for the presented calculation [where λ=λ(Q) and r≠r(Q)]. The use of separate equations for hydraulic resistance, r or at least separate comments for the laminar and turbulent regimes in Eqs. (1-3) of the discussed paper should be avoided (better to say the parameter, λ should always be used and not r). However, in laminar, turbulent zone and also for the unstable transient zone between them (Table 2 and4 Finally, the flow dependent Jacobian, as recommended by [START_REF] Simpson | Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model[END_REF], and the Darcy-Weisbach scheme should be used. Furthermore, EPANET also should be revised in light of the following discussion. In the laminar zone the Darcy friction factor should be calculated using λ=64/Re. For the turbulent zone, the Swamee and Jain approximation [used in EPANET, in the discussed paper and in [START_REF] Spiliotis | Water distribution system analysis: Newton-Raphson method revisited[END_REF]] should be replaced with one of the more accurate approximations of the Colebrook equation [where λ=λ(Re, ε/D)]. The Swamee and Jain approximation produces maximal error of 2% compared to the iterative solution of the Colebrook equation as shown in the recent study by Brkić (2011a).

For example, the approximation proposed by [START_REF] Romeo | Improved explicit equation for estimation of the friction factor in rough and smooth pipes[END_REF] with a maximal error up to 0.1345% is more accurate. Also, approximations by [START_REF] Buzzelli | Calculating friction in one step[END_REF] and [START_REF] Serghides | Estimate friction factor accurately[END_REF] or by [START_REF] Zigrang | Explicit approximations to the solution of Colebrook friction factor equation[END_REF] are very accurate with an error of no more than 0.1385%. [START_REF] Vatankhah | Discussion of "Turbulent flow friction factor calculation using a mathematically exact alternative to the Colebrook-White equation[END_REF] also show very accurate approximation with a relative error up to 0.1472%. There are at least ten approximations more accurate than the Swamee-Jain equation. These very accurate approximations are sometimes complex but they can be simply included in a computer code (Matlab, Excel, etc.) without causing significant additional computational time. Also, between the laminar and turbulent zones it is more convenient to model it as a sharp jump (sudden transition to turbulence) instead of using a cubic spline (explained in the discussed paper in Table 2; also instead of a cubic spline, the simple equation λ=0.0025•Re 0.33 is more convenient). At last but not less important, [START_REF] Churchill | Friction-factor equation spans all fluid flow regimes[END_REF] also provides a single formula which covers all three cases in the discussed paper (Table 2 in the discussed paper of [START_REF] Simpson | Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model[END_REF], i.e. used for the laminar and turbulent zone as well as for the unstable transient zone between them as described in EPANET). The equation by [START_REF] Churchill | Friction-factor equation spans all fluid flow regimes[END_REF] perfectly fits the laminar regime and provides a reasonable transition to turbulence (cubic spline in the discussed paper). Finally it approximates the implicit Colebrook function (suitable only for the turbulent zone) at least equally as accurate as the Swamee-Jain approximation. [START_REF] Simpson | Jacobian matrix for solving water distribution system equations with the Darcy-Weisbach head-loss model[END_REF] use Matlab and hence for the readers" convenience, Matlab code for the well known equation by [START_REF] Churchill | Friction-factor equation spans all fluid flow regimes[END_REF] is given:

  approximations of Colebrook"s equation shown inBrkić (2011a) are more suitable for the calculation of the derivatives used in the Jacobian matrix since the original Colebrook equation is implicit in the Darcy friction factor, λ. The implicitly given form of the original Colebrook function complicates the calculation of the first derivative or even makes it impossible. The Reynolds number in the flow term Re=4•Q/(D•π•ν) has to be placed in the equation for the Darcy friction factor, λ (the Swamee and Jain approximation in Eq. 17) and the first derivative where flow, Q, being variable, is treated as the variable and calculated from Eq. (17).

  of the discussed paper) it is certain that the Darcy friction factor, λ depends on the Reynolds number [λ=λ(Re, ε/D)] and subsequently in general on flow [where Re=Re(Q)]. This is not true for the fully turbulent regime (fully developed flow) where the Darcy friction factor, λ depends only on relative roughness ε/D and not on the Reynolds number Re and where, therefore the influence of the flow Q can be neglected).

  function Lambda = churchilllamapp(Re,RelEpsilon) Lambda = (2.457.*ln(1./((7./Re).^0.9+0.27.*RelEpsilon))).^16; Lambda = 8.*((8./Re).^12+(1./(Lambda+(37530./Re).^16)^1.5)).^(1/12); Clamond (2009) provides Matlab code for the iterative solution of the implicit Colebrook equation. Acknowledgment I appreciate the assistance of prof. Žarko Ćojbašić from the University of Niš (Serbia) who provides the valuable help with the Matlab code. Sincere thanks Mr. Simpson for his prompt response to my e-mail inquiry about here discussed paper. Notation The following symbols are used in this discussion: D -inner pipe diameter (m) Qflow rate (m 3 /s) Δppressure drop (Pa) Lpipe length (m) λ -Darcy (Darcy-Weisbach) or Moody friction factor (-) ["Lambda" in the presented Matlab code] rhydraulic resistance defined in text ε/Drelative roughness of inner pipe surface (-) ["RelEpsilon" in the presented Matlab code] Re -Reynolds number (-) ["Re" in the presented Matlab code] νkinematic viscosity (m 2 /s) ρdensity (kg/m 3 ) π ~ 3.1415


for presenting a possible way to establish sustainable electric energy system in Croatia.

The discusser, however, would like to add few points and some literature references that might be of interest for readers.

. As seen in Boţiĉević Vrhovĉak et al. (2006), total installed capacity in Croatian electric power system is 4GW (52% in hydro, 29% oil, 8% in nuclear, 8% in coal and 3% in gas). Electricity consumption is about 15 TWh (Schneider et al., 2007). Croatia has a strong but not very efficient co-generation sector,

delivering 12% of the final energy [START_REF] Lipošćak | Sustainability assessment of cogeneration sector development in Croatia[END_REF][START_REF] Lonĉar | An analysis of the legal and market framework for the cogeneration sector in Croatia[END_REF][START_REF] Uran | Feed-in tariff and market electricity price comparison: The case of cogeneration units in Croatia[END_REF]. As in all transition countries, electricity consumption in Croatia sharply declined after 1991. The reason for this in Croatia was not only the breakdown of the former communist regime but also civil war in former Yugoslavia that lasted from 1991 to 1995. For a decade now, Croatia has been experiencing increased electricity consumption, estimated to be growing by approximately 500GWh annually [START_REF] Schneider | Mapping the potential for decentralized energy generation based on renewable energy sources in the Republic of Croatia[END_REF]. The net electricity import is 31%-40% of the total consumption [START_REF] Schneider | Mapping the potential for decentralized energy generation based on renewable energy sources in the Republic of Croatia[END_REF][START_REF] Saner | Is it only CO 2 that matters? A life cycle perspective on shallow geothermal systems[END_REF]). About 0.421kg/kWh of CO 2 is released in atmosphere by electricity production in Croatia, where this value for continental Europe is estimated to 0.599kg/kWh of CO 2 [START_REF] Saner | Is it only CO 2 that matters? A life cycle perspective on shallow geothermal systems[END_REF]. Compared to the energy consumption collapse in some transitional countries, Croatia has passed through a relatively short-term reduction of GHG emissions since 1990 because of higher energy efficiency of its pre-transition economy [START_REF] Duić | Croatia energy planning and Kyoto Protocol[END_REF]. The external costs caused by human health impacts of power plants in Croatia are comparable to those attached to the EU power plants (Boţiĉević Vrhovĉak et al., 2005).

All regions in Croatia are not equally suitable for certain renewable energy source. Significant differences exist between northern region (with part of Pannonian plain) and southern Adriatic region. Total coastal length of Croatia, including over 800 islands and islets, equals 5790 km.

Balancing tool for islands that simulates the integration of renewable energy into energy-systems H 2 RES was developed in 2000 by the Instituto Superior Técnico in Lisbon, Portugal and the Faculty of Mechanical Engineering and Naval Architecture in Zagreb, Croatia [START_REF] Duić | RenewIslands methodology for sustainable energy and resource planning for islands[END_REF][START_REF] Krajaĉić | H 2 RES, Energy planning tool for island energy systems-The case of the Island of Mljet[END_REF]). H 2 RES balances the hourly time series of water, electricity, heat, and hydrogen demand, as well as appropriate storages and supply over any user-defined period. The tool has been specifically designed to increase the integration of renewable sources and hydrogen, into island energy-systems which operate as stand-alone systems. It can serve as a planning tool for single wind, hydro, or solar power systems, as well as for planning larger energy-systems.

The current share of renewable energy sources in electricity production in Croatia is high, around 50% (Boţiĉević Vrhovĉak et al., 2006). The potential of renewable energy sources for electricity generation still has to be carefully analyzed. The estimated potential varies between 3.5TWh in a moderate projection to 6TWh per year in the most optimistic perspective (Boţiĉević Vrhovĉak et al., 2006). In the year 2000, total primary energy consumed was about 360PJ [START_REF] Karasalihović | Natural gas in Croatia's energyfuture[END_REF], out of which 75PJ was produced from renewable sources (Boţiĉević Vrhovĉak et al., 2006). In 1997, the Croatian government started five national energy programs dealing with renewable energy sources: BIOEN (biomass and waste), SUNEN (solar energy), ENWIND (wind energy), MAHE (small hydro) and GEOEN (geothermal energy). Within the national energy program BIOEN, it has been demonstrated that by 2030 at least 15% of the total energy consumed could be obtained from biomass [START_REF] Schneider | Mapping the potential for decentralized energy generation based on renewable energy sources in the Republic of Croatia[END_REF]. More than 40% of Croatian territory is covered by forests, which means that the potential for biomass use is significant (Ćosić et al., 2011). Small hydro power plants (<10MW) can provide 2.37GWh annually of electrical energy in Croatia [START_REF] Schneider | Mapping the potential for decentralized energy generation based on renewable energy sources in the Republic of Croatia[END_REF]. The use of wind energy in Croatia is currently very low, but future wind generators at the best sites in Croatia are close to becoming marginally competitive with fossil-fuel technologies [START_REF] Feretić | Feasibility analysis of wind-energy utilization in Croatia[END_REF][START_REF] Feretić | Probabilistic analysis of electrical energy costs comparing: production costs for gas, coal and nuclear power plants[END_REF]. On the Adriatic basin, in some area along the coast, a wind power of over 300W/m 2 and more than 2500kWh/m 2 per year may be gained with potential estimated at 200-400MW installed capacity and about 800GWh yearly production [START_REF] Poje | Assessment of wind energy potential in Croatia[END_REF][START_REF] Urli | Hybrid photovoltaic/wind grid-connected power plants in Croatian renewable energy program[END_REF][START_REF] Schneider | Mapping the potential for decentralized energy generation based on renewable energy sources in the Republic of Croatia[END_REF]. Solar potentials in Croatia with 2300-2800 sunshine hours per year in particular Adriatic islands can be estimated up to 5.5GJ/m 2 per annum with peak at Dubrovnik area (Desnica et al. 1988[START_REF] Sinjeri | Solar radiation on variously oriented collectors in Croatia[END_REF][START_REF] Urli | Hybrid photovoltaic/wind grid-connected power plants in Croatian renewable energy program[END_REF][START_REF] Franković | Thermal utilisation of solar energy: Croatian national energy program SUNEN, 1998-2010[END_REF][START_REF] Hrastnik | Solar energy demonstration zones in the Dalmatian region[END_REF]. The potential of already discovered geothermal fields is 839MW of thermal energy and 30-47.9MW of electric energy [START_REF] Schneider | Mapping the potential for decentralized energy generation based on renewable energy sources in the Republic of Croatia[END_REF][START_REF] Guzović | Possibilities of electricity generation in the Republic of Croatia by means of geothermal energy[END_REF]. The results of the computational study by [START_REF] Davis | Geothermal power production from abandoned oil wells[END_REF] showed that abandoned oil wells (also exist in Croatia) have the potential to produce a significant amount of power when they are modified to become double-pipe heat exchangers. The calculations [START_REF] Davis | Geothermal power production from abandoned oil wells[END_REF] show that the net power produced from such a well may exceed 3MW for a bottom-hole temperature of 450K and an injection pressure of 30x10 5 Pa. This amount of power is not intermittent as with other renewable energy sources, and it may be available for peak demand as well as basic demand. The power produced varies and depends significantly on the down-hole temperature, the injection pressure, the injection velocity, and the geometric characteristics of the pipe (inner and outer radius and insulation thickness).

Electricity for hot water production with share of 43% represents significant share of the total consumption in Croatia (Atanasiu and Bertoldi, 2008). Electricity is used also for heating and this amount is significant in the residential consumption. Electricity for heating can be substitute with natural gas (Brkić and Tanasković, 2008). This can be done especially in northern part of Croatia, where natural gas supply system has been sufficiently developed. Natural gas can be significant if Croatia became transient country for natural gas from Russia to market of EU [START_REF] Brkić | Serbian gas sector in the spotlight of oil and gas agreement with Russia[END_REF]. Natural gas can be also implemented at the Adriatic coast and islands [START_REF] Vujĉić | The status and prospects for the gasification in the County of Split and Dalmatia[END_REF][START_REF] Karasalihović | Natural gas in Croatia's energyfuture[END_REF][START_REF] Dondi | Croatia methane pipeline project requires difficult canal crossing[END_REF]. Energy efficiency is also very important for control of electric power consumption. Gas domestic boilers and electrical appliances with higher class of efficiency should be used (Brkić and Tanasković, 2008). Of course, district heating from remote plant in towns fired by gas can be also implemented. Share of sales of A class appliances, including domestic gas boilers in Croatia is 26.6% (Atanasiu and Bertoldi, 2008). The sales of the main domestic appliances in Croatia have a still higher share of lower energy classes. In Croatia was registered the lowest share of A class washing machines with 17.6% sold on the market (and only 0.8 of A+ class) comparing with the average share of 71.2% A and 8.4% A+ class from the EU markets (Atanasiu and Bertoldi, 2008). Building materials (e.g. windows) have become cheaper over the last decades in spite of significantly increasing energetic quality. Now, costs for building materials are 12% lower in Croatia compared with Germany.

References:

Atanasiu, B., and Bertoldi, P. (2008). "Residential electricity consumption in New Member States and Candidate Countries." Energ. Buildings, 40(2), 112-125. Boţiĉević Vrhovĉak, M., Tomšić Ţ., and Debrecin N. (2005). "External costs of electricity production: case study Croatia." Energ. Policy, 33(11), 1385-1395. Boţiĉević Vrhovĉak, M., Tomšić Ţ., and Debrecin N. (2006). "Potential and use of renewable energy sources in Croatia." Renew. Energ., 31(12), 1867-1872. Brkić, D., and Tanasković, T.I. (2008). "Systematic approach to natural gas usage for domestic heating in urban areas. " Energy, 33(12), 1738" Energy, 33(12), -1753. .


Dejan Brkić
email: dejanrgf@tesla.rcub.bg.ac.rs

Vladimir Mitrović

. The discussers are also involved in teaching of underground flow to the students of petroleum engineering, but they use different approach; i.e. students have to measure some quantities in the laboratory and then to calculate others.

Petroleum department at University of Belgrade had very simple professional permeameter, but few parts are now broken and since the spare parts are not available, the discussers have to manage to made exercises with student using hand-made apparatus. Experience shows that student cannot understand concept of hydraulic conductivity and permeability before few introductive examples. These exercises will be used also to illustrate problem which can arise

with use of corn syrup as working fluid. Simple system of sphere packing is used to calculate porosity, permeability and specific surface of porous systems (Figure 1). Most models of flow through porous media are first constructed with perfect spheres in mind. Note that system is spatial (where granular beads are spheres and not circles). Also, different colors of spheres used in Figure 1 do not necessary mean that particles have different mineralogical structure. Color is used only to differentiate particles by their size (or precisely by diameter). reservoirs where pressure overburden is significant [START_REF] Trussell | Review of flow through porous media as applied to head loss in water filters[END_REF]. First task for student is to calculate porosity using knowledge from geometry (answers: Φ A) =Φ B) ≈0.477 and Φ C) =Φ D) ≈0.26). For Darcy's bottle more possible pattern is A), i.e. B). Example E) is introduced to examine mixture of small and large beads packed similar as in example A) or B). It is not always obvious at first sight for students that Φ A) =Φ B) , but Φ E) <Φ A) and consequently Φ E) <Φ B) .

Next task is to calculate permeability. This can be done using Kozeny-Carman model [START_REF] Reyssat | Imbibition in layered systems of packed beads[END_REF] for spherical systems (1):

  23 2 180 1 d k    (1)
Fixing value of porosity, it is obvious that larger diameter of beads produces increased permeability (k A) >k B) and k C) >k D) ). Third task is to calculate specific surface which is for a media made of uniform spheres well defined (2):

  61 s d     (2)
Here, specific surface is defined as the surface area of the media divided by its bulk volume where s A) <s B) and s C) <s D) . Definitions of other types of specific surfaces are available from paper of [START_REF] Trussell | Review of flow through porous media as applied to head loss in water filters[END_REF]. Porosity, permeability and specific surface are defined only using syrup is not adequate fluid because then concept of wettability (adsorption effect) has to be introduced in this case. The surface wettability of porous system is the result of complex interactions between fluids and the surface of sand or gravel grains. Measurements of contact angles (Figure 2) provide a convenient and conceptually simple approach to quantification of the wettability given by liquid pairs at a smooth mineral surface.

Figure 2. Contact angles caused by different wettability cause changes in effective porosity

It is clear that water and corn syrup has different wetability capacities. Can one prove that permeability does not depend on fluid properties only knowing hydraulic conductivity for water and corn syrup (can be measured using Darcy's bottle), and knowing water and corn syrup kinematic viscosity (3):

w CS air w CS air k K K K g g g                            (3)
Previous equation is valid only in absence of wetting effect (thus can be accomplished using air as working fluid or in theory even mercury). According to Benner and Bartell (1941), even water on same mineral surfaces can produce different wetting angles in presence of different fluids (Figure 3). Mineral surface in examples 1) and 2) is quartz and in 3) and 4) is calcite. Water is surrounded with i-octane in example 1) and 3) and with naphthenic acid in example 2) and 4). 4). Residual of corn syrup in sample is greater than water residual. This means that increased amount of non-movable corn syrup in comparison with non-movable water cause effect of decreased effective porosity and consequently not only decreased value of Darcy's conductivity but equally of permeability (similar e.g. one can conclude from eq. y=A/B=3/4 that A=3 and B=4, but if A=sin(x) and B=cos(x) previous conclusion fails). This amount of nonmovable water can be extracted entirety from the sand sample by drying, but with corn syrup this cannot be done. For example, air permeability in presence of water varied from 8.5·10 -12 m 2 for volumetric water content of 2.03% to 2.3·10 -12 m 2 for water content of 20.0% [START_REF] Springer | Air permeability of porous materials under controlled laboratory conditions[END_REF]. It is important to note that Darcy's experiments actually provide no information about any properties within a packed column [START_REF] Gray | Examination of Darcy's law for flow in porous media with variable porosity[END_REF]. For measure of Darcy's conductivity for teaching purpose, we can recommend apparatus (Figure 5) developed at the University of Washington [START_REF] Massmann | Exercises illustrating flow in porous media[END_REF]. Using air as working fluid, sand remains dry (but for very accurate results Klinkenberg's effect has to be resolved). For minor pressure drop air can be treated as incompressible fluid.

Figure 5. Simple hand-made permeameter

The conductivity of the sample (Figure 5) is measured by submerging the lower section of the tube in water and measuring how quickly the water forces the air through the soil sample. First, the top of the tube is sealed and then cork is removed. Because the air pressure in the lower section of the tube is greater than atmospheric pressure water level in the tube will arise. Time in which water level is changed from H 0 to H 1 also has to be measured (4).

            (4)
Further, permeability and Darcy's conductivity for water and corn syrup can be calculated using Eq. (3). [START_REF] Carlton | Home-made head tanks can be used to demonstrate Darcy's law[END_REF] also proposed simple way to make Darcy's bottle. Permeameter can be made also using instructions in the paper of [START_REF] Noll | Building bridges between field and laboratory studies in an undergraduate ground water course[END_REF].

Figure 1 .

 1 Figure 1. System of sphere packing used to calculate porosity, permeability and specific surface

  properties of sample matrix, but the value of hydraulic conductivity in Darcy's equation (Eq. 1 of the original paper) incorporates the characteristics of the porous media as well as the fluid. It has unit of length/time (m/s). In correcting Darcy's coefficient for the viscosity and the density of the fluid, permeability eliminates the influence of the type of fluid on the permeability which are constant for the model, making the model more robust. Permeability has units of lengthsquared (m 2 ). Correlation of Darcy's hydraulic conductivity and permeability can be calculated using kinematic viscosity of fluid as k=K·(ν/g). Here lies first problem in experiment proposed with Darcy's bottle. Darcy's conducted experiments with water and not with other type of fluids. Corn

Figure 3 .

 3 Figure 3. Different water wetting angles on same mineral surfaces in presence of different fluids

Figure 4 .

 4 Figure 4. Process of filling and empting a Darcy's bottle

  .C., and Bartell, F.E. (1941). "The effect of polar impurities upon capillary and surface phenomena in petroleum production." Drill.Prod. Prac. API,[341][342][343][344][345][346][347][348]

Figure 1 .Figure

 1 Figure 1 DB Click here to download high resolution image


Dejan Brkiđ
email: dejan.brkic@jrc.ec.europa.eu

Keywords: Colebrook equation, ColebrookͲWhite friction factor, Hydraulic resistance, Error analysis, Turbulent flow, Pipes

Maximal relative error of the explicit approximation to the Colebrook equation for flow friction presented in the discussed paper by Saeed Samadianfard [J. Pet. Sci. Eng. 92Ͳ93 (2012), 48Ͳ55; doi. 10.1016/j.petrol.2012.06.005] is investigated. Samadianfard claims that his approximation is very accurate with the maximal relative error of no more than 0.08152%. Here is shown that this error is about 7%. Related comments about the paper are also enclosed.

1. Introduction Samadianfard (2012) presents an approximate explicit formula as a replacement for the implicitly

H O O

(1) Equation (1) [Eq. 1a of the discussed paper] is from [START_REF] Colebrook | Turbulent flow in pipes with particular reference to the transition region between the smooth and rough pipe laws[END_REF] but where coefficient 3.71 is slightly different, i.e. 3.7. This unknown friction factor, ʄ 0 , cannot be extracted to be on the left side of the equal sign analytically, i.e. with no use of some simplifications. Better to say, it can be expressed explicitly only if some kind of approximate calculus takes place, where ʄ 0 o ʄ. Samadianfard (2012) (2) Samadianfard (2012) used genetic programming to develop his approximation. ojbašiđ and Brkiđ (2013) also used genetic algorithms to improved accuracy of some of the available approximations [START_REF] Romeo | Improved explicit equation for estimation of the friction factor in rough and smooth pipes[END_REF][START_REF] Serghides | Estimate friction factor accurately[END_REF].

2. Error analysis Samadianfard (2012) claims that his approximation (2) [Eq. ( 29) of the discussed paper] produces maximal relative error, ɷ max of no more than 0.08152% compared with the very accurate iterative solution of the implicit Colebrook equation, ʄ 0 . Using the procedure from Brkiđ (2011a), it can be shown that this error, ɷ max in the practical range of applicability of the Reynolds number (Re) and the relative roughness of inner pipe surface (ɸ/D) is about 7.4289%1 [mean (average) relative error is about 1.6198%]. Distribution of the relative error of approximation by Samadianfard (2012) [(2); Eq.

(29) of the discussed paper] can be seen in Figure 1. Relative error is defined as ɷ max = (ʄ 0 Ͳ ʄ)/ʄ 0 ͼ100% as shown by Eq. ( 26) of the discussed paper where ʄ 0 is very accurate iterative solution of the implicit Colebrook equation while ʄ is the solution of the approximation by Samadianfard ( 2012) [(2); Eq. ( 29) of the discussed paper]. The error is evaluated twice in MS Excel according to methodology from Brkiđ (2011a[START_REF] Brkiđ | Determining friction factors in turbulent pipe flow[END_REF]. To produce Figure 1, entire practical domain of the Reynolds number (Re) and the relative roughness of inner pipe surface (ɸ/D) is covered with 740 pointͲmesh. Using double check in MS Excel, second mesh of 860 points is formed using 20 values of the relative roughness (ɸ/D) shown in Table 1 [Table 2 of discussed paper] and using 43 values of the Reynolds number (Re); from 4ͼ10 3 to 10 4 with pace 10 3 , from 10 4 to 10 5 with pace 10 4 , from 10 5 to 10 6 with pace 10 5 , from 10 6 to 10 7 with pace 10 6 , and from 10 7 to 10 8 with pace 10 7 . Similarly as in Samadianfard (2012), results are shown in Table 1 and compared with the results from Table 2 of the original paper by Samadianfard (2012). According to this second check, the maximal relative error ɷ max is about 6.9107% and mean (average) relative error in the range of applicability of equation is about 1.2595%. Mean square error as defined by Eq. 27 of the discussed paper, for the mesh of 740 points is about 1.1460ͼ10 Ͳ7 and for the mesh of 860 points is about 2.2711ͼ10 Ͳ7 . According to [START_REF] Winning | Explicit friction factor accuracy and computational efficiency for turbulent flow in pipes[END_REF], the approximation by Samadianfard (2012) is in the group with medium value of mean square error (very small is lower than 10 Ͳ11 , small is between 10 Ͳ11 and 10 Ͳ8 , medium is between 10 Ͳ8 and 5ͼ10 Ͳ6 , and large is above 5ͼ10 Ͳ6 ). Samadianfard (2012) reported value of mean square error of 1.95ͼ10 Ͳ10 .

Source of the error

Regarding the error from Table 1 and from Figure 1 of this discussion, it is useful to mention that actually three level of the accuracy can be introduced; 1. Colebrook's equation is empirical (other maybe more accurate equations can be used to describe related physical processes), 2. Colebrook's equation can be solved very accurately using iterative procedure (term "accurate by default" can be used or better to say, this error can be neglected in many cases; this is ʄ 0 in this discussion), 3.

Relevant explicit approximations of the Colebrook equation (this is ʄ in this discussion) can be used to avoid iterative procedure (their error can be estimated very accurately compared with the error of iterative solution). The estimated error from Table 1 and from Figure 1 of this discussion belongs to the category explained in point 3 (where friction factor from point 2 assumed as accurate). Using "philosophy" from point 2, value ʄ 0 of friction factor is calculated, while friction factor ʄ is calculated using point 3. Finally, the maximal relative error was calculate as ɷ max = (ʄ 0 Ͳʄ)/ʄ 0 ͼ100%. It is also useful to note that the Moody diagram [START_REF] Moody | Friction factors for pipe flow[END_REF]) is only graphical interpretation of the Colebrook equation [START_REF] Colebrook | Turbulent flow in pipes with particular reference to the transition region between the smooth and rough pipe laws[END_REF]). Moody's diagram is not less or more accurate compared with Colebrook's equation. So, inaccurate reading from Moody's diagram can introduce one more source of error [START_REF] Yŧldŧrŧm | ComputerͲbased analysis of explicit approximations to the implicit ColebrookͲ White equation in turbulent flow friction factor calculation[END_REF][START_REF] Brkiđ | New explicit correlations for turbulent flow friction factor[END_REF][START_REF] Fang | New correlations of singleͲphase friction factor for turbulent pipe flow and evaluation of existing singleͲphase friction factor correlations[END_REF], which is maybe possible cause of the value of the error estimated in Samadianfard (2012). In this case, source of the error is not then calculative, but rather systematic, where the main source of the error is in inaccuracy of ʄ 0 (which should be accurate by default).

Analysis of the equation

It is worth to mention that the approximation by Samadianfard (2012) does not contained a single logarithmic expression. This structure of the equation is valuable contribution since the logarithmic expression uses a lot of computational resources in a computer environment [START_REF] Giustolisi | Some explicit formulations of Colebrook-White friction factor considering accuracy vs. computational speed[END_REF][START_REF] Clamond | Efficient resolution of the Colebrook equation[END_REF]). Only, three other available approximations do not use logarithms; [START_REF] Moody | An approximate formula for pipe friction factors[END_REF], [START_REF] Wood | An explicit friction factor relationship[END_REF] and [START_REF] Chen | A simple explicit formula for the estimation of pipe friction factor[END_REF], and all three are less accurate than the approximation by Samadianfard (2012). But also, Samadianfard (2012) uses nonͲinteger power such as Re ɸ/D which in computer environment usually means exp[(ɸ/D)ͼln(Re)] where 'ln' is natural (Napierian) logarithm [START_REF] Clamond | Efficient resolution of the Colebrook equation[END_REF].

Reference notes

Finally, it has to be mentioned that source of the review article of Brkiđ (2011a) is not Petroleum Science and Technology [ISSN 1091Ͳ6466 (Print), 1532Ͳ2459 (Online)] as Samadianfard (2012) Also, it can be mentioned that Colebrook equation is published solely by [START_REF] Colebrook | Turbulent flow in pipes with particular reference to the transition region between the smooth and rough pipe laws[END_REF] and not in collaboration with White as often cited. [START_REF] Colebrook | Turbulent flow in pipes with particular reference to the transition region between the smooth and rough pipe laws[END_REF] acknowledged contribution of White in a footnote of his paper with the reference of their previous joint work [START_REF] Colebrook | Experiments with fluid friction in roughened pipes[END_REF].

Notation remark

In this discussion, the Darcy friction factor is labeled as ʄ while in the discussed paper as f. This is because f should be rather used for Fanning friction factor. This can make confusion since the physical meaning is equal but where the Darcy friction factor is 4 times greater than Fanning's friction factor.

Conclusion

According to the classification proposed by Samadianfard (2012), his approximation with the relative error of about 7% (7.43% after first evaluation and 6.91% after second evaluation; both in MS Excel), have to be characterized as less accurate method (group of methods with relative error of more than 3%) and should be placed in Figure 3 of the discussed paper (or even in Figure 2 of the discussed paper with the nonͲadvisable approximations) and not in Figure 6 of the discussed paper with the extremely accurate approximations.

given

  Colebrook equation for fluid flow pipe friction. The empirical Colebrook equation (1) [Eq. 1a of the discussed paper] relates the unknown flow friction factor (ʄ 0 ) with the known Reynolds number (Re) and the known relative roughness of inner pipe surface (ɸ/D). Widely used empirical and nonlinear Colebrook's equation for calculation of Darcy's friction factor is iterative i.e. implicit in fluid flow friction factor since the unknown friction factor appears on the both sides of the equation [ʄ 0 =f(ʄ 0 , Re, ɸ/D)] (Colebrook 1939).

  presents one explicit approximation of the Colebrook equation [ʄуf(Re, ɸ/D)]. His explicit formula is given as (2) [Eq. (29) of the discussed paper].

Figure 1 .

 1 Figure 1. Distribution of the relative error of approximation by Samadianfard (2012) [(2); Eq. (29) of

  cited in the discussed paper but Journal of Petroleum Science and Engineering [ISSN 0920Ͳ4105]. In the first mentioned journal Brkiđ (2011c) published an article with his explicit approximation of Colebrook's equation, which possible makes confusion here.

Table 1 .

 1 Distribution of maximal relative error of the approximation by Samadianfard (2012) [(2); Eq.

	(29) of the discussed paper]

Table 1 .

 1 Distribution of maximal relative error of the approximation by Samadianfard (2012) [(2); Eq. (29) of the discussed paper]

		Relative error (%)
		Samadianfard	This discussion
	Relative roughness (ɸ/D)	Max ɷ max	Mean ɷ	Max ɷ max	Mean ɷ
	1ͼ10 Ͳ6	0.049	0.012	6.3210	2.9289
	5ͼ10 Ͳ6	0.063	0.02	6.9107	3.0483
	1ͼ10 Ͳ5	0.07	0.024	6.2389	3.1781
	5ͼ10 Ͳ5	0.076	0.035	4.6612	2.7105
	1ͼ10 Ͳ4	0.077	0.039	3.9580	2.0012
	2ͼ10 Ͳ4	0.078	0.043	3.1845	1.3767
	4ͼ10 Ͳ4	0.078	0.047	2.1375	0.9998
	6ͼ10 Ͳ4	0.078	0.049	1.3145	0.7962
	8ͼ10 Ͳ4	0.079	0.051	1.1123	0.6961
	1ͼ10 Ͳ3	0.079	0.052	1.2409	0.6857
	2ͼ10 Ͳ3	0.079	0.055	2.4851	1.1169
	4ͼ10 Ͳ3	0.079	0.06	4.0676	1.1633
	6ͼ10 Ͳ3	0.079	0.062	4.6207	1.0094
	8ͼ10 Ͳ3	0.079	0.065	4.6890	0.8434
	1ͼ10 Ͳ2	0.079	0.066	4.5401	0.6976
	1.5ͼ10 Ͳ2	0.079	0.07	3.8766	0.5606
	2ͼ10 Ͳ2	0.079	0.072	3.1906	0.4685
	3ͼ10 Ͳ2	0.079	0.074	2.0914	0.3450
	4ͼ10 Ͳ2	0.079	0.075	1.2864	0.2831
	5ͼ10 Ͳ2	0.08	0.076	0.6478	0.2804

Graphical Abstract (for review) DB

Figure

1 DB

This discussion was initially submitted on 8 December

2013, and on 12 December 2013 the discussion of Vatankhah (2014) appeared onͲline with similar findings regarding this error analysis