Mathematical model on the GasModels.jl manual

[cvr]

Hi, One question regarding the math in the online documentation. Before integration, when writing p dp/dx = - λ a² ϕ |ϕ| / (2 D) it should be equal to writing d(p²)/dx = - λ a² ϕ |ϕ| / D. Thus, after integration you should get p²(L) - p²(0) = - L λ a² ϕ |ϕ| / D = - L λ a² f |f| / (A² D), however you have an extra ½ in your equations. Am I missing something? If not, the implementation in the code is affected by this ½ factor or is just a glitch in the manual.

Thank you for your work on GasModels.jl and your decision on releasing the code. I'm planning on using it for lecturing purposes, if there is no conflict.

Best regards.

[kaarthiksundar]

@cvr Sorry for the confusion, the derivation of equation is not well explained in the online documentation. Please take a look at our paper https://arxiv.org/pdf/1803.07156.pdf for a detailed derivation (see page 5). The confusion arises from the RHS of the equation p²(L) - p²(0) = - L λ a² ϕ |ϕ| / D. Here ϕ |ϕ| is actually average flow in the pipe i.e. (ϕ(0)+ ϕ (L)) |ϕ(0)+ ϕ(L)|. Then your derivation will go through. If it is still not clear, let us know and we shall update the documentation to reflect the paper.

[cvr]

Thank you for your reply and the paper.

[ccoffrin]

@kaarthiksundar, maybe we should add a reference to the arxiv paper to the docs and then we can close this issue?

[rb004f]

After some discussion, we believe the suggested change is a better approximation in a purely steady-state regime with long pipe distances. We will adjust the formulation shortly,

[cvr]

I have analysed the equations and I believe the ½ factor in the equations is not needed, as well as the need of computing the average mass flow as stated by @kaarthiksundar. Sorry for the lengthy discussion that follows, but it's difficult to be synthetic.

Assuming two nodes (several mass flows entering/exiting) and one pipe between them, for the pipe itself there is only one inlet and outlet. Lets further assume the pipe's diameter is not constant between inlet and outlet. Applying the mass conservation equation in its full 3D form, ∂ρ/∂t + ∇·vec(ϕ) = 0, and integrating it over the volume yields ∂/∂t(∫∫∫ρ dV) + Σvec(ϕ)·vec(A) = 0. If there are no mass sink/sources throughout the pipe (e.g. leakage), then even for compressible flows with sub-to-supersonic transitions the mass gain/loss only happens at the inlet/outlet, thus: Σvec(ϕ)·vec(A) = ΣϕA = Σf = 0 which for the pipe means f_in = f_out. As such, the average mass flow through the pipe must be always equal to the mass flow at the inlet and at the outlet. The average may be performed but it should yield f_med = (f_out + f_in)/2 = 2 f_in / 2 = f_in.

If you refer to Chapter 2 in "Evaluating Gas Network Capacities", equation (2.26-2.17) (and the more generic form (2.24)) shows p²(L) - p²(0) = - L λ a² f |f| / (A² D), so without the ½ factor.

I presume some mix-up may arise if the continuity eq. is given purely in 1D form, ∂ρ/∂t + ∂ϕ/∂x = 0, overlooking its original 3D form. Even for ρ=cte, this simplified assumption is only valid for constant diameter pipes, as direct integration yields ϕin = ϕout => fin = fout Ain / Aout, yielding different mass flux at the inlet and outlet (which cannot be). The misjudgement is that for expanding or contracting pipes the flow streamlines are straight. There are streamlines in the y and z directions according to the area changes. While the 3D form of the equations consider these (even if the end result simplifies to 1D), a pure 1D form of the continuity equation will fail to represent this.

In engineering fluid mechanics, the manometric head loss due to friction is represented by the Darcy-Weisbach equation as Hfriction = U² λ L / (2 g D) and local losses through Hlocal = U² K / (2 g). These are related to the fluid velocity through the loss, which changes with area and density. For a 1D pipe composed of several sections (areas or diameters) and local losses (valves, curves, etc), one usually substitutes U = f / (ρ A). As the mass flux is constant (independent if the fluid is compressible or not) one may get Hlosses = f |f| / (2 g) [Σi fi Li/(Di ρi² Ai²) + Σj Kj/(ρj² Aj²)]. This is a condensed way to get all losses between a system with a single inlet and outlet.

Best regards

[rb004f]

Thank you @cvr . We are in agreement and we have adjusted the model accordingly.